AD-A240  656 


NAVAL  POSTGRADUATE  SCHOOL 


Monterey,  California 


THESIS 


-4 


|e:lecte 

S£a23..1991] 


DISCRETE  ARMA  MODEL  FOR  NATURAL 
RESONANCES  IN  ELECTROMAGNETIC  AND 
ACOUSTIC  SCATTERING 

by 

Yuval  Cohen 
September  1990 

Thesis  Advisor:  Michael  A.  Morgan 


Approved  for  public  release;  distribution  is  unlimited. 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


la  REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


2a  SECURITY  CLASSIFICATION  AUTHORITY 


2b.  DECLASSIFICATION /DOWNGRADING  SCHEDULE 


A  PERFORMING  ORGANIZATION  REPORT  NUMBERtS/ 


REPORT  DOCUMENTATION  PAGE 


lb  RESTRICTIVE  MARKINGS 


Form  Approved 
0MB  No.  0704-0188 


3  DISTRIBUTION /AVAILABILITY  OF  REPORT 

Approved  for  public  release; 
distribution  is  unlimited 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


6a  NAME  OF  PERFORMING  ORGANIZATION  6b  OFFICI 

(If  ap 

Naval  Postgraduate  School  EC 


6c.  ADDRESS  {City,  State,  and  ZIP  Code) 

Monterey,  CA  93943-5000 


6b  OFFICE  SYMBOL  7a..NAME  OF  MONITORING  ORGANIZATION  " 
(If  applicable) 

EC  Naval  Postgraduate  School 


7b.  ADDRESS  (C/fy,  State,  and  ZIP  Cod&) 

Monterey,  CA  93943-5000 


8a  NAME  OF  FUNDING  .  SPONSORiNC 
ORGANIZATION 


8c.  ADDRESS  (City,  State,  and  ZIP  Code) 


8b.  OFFICE  SYM80L  9  PRGCUSEMENT  INSTRUMEN-T^-IDENTIFICATION  NUMBER 
(If  applicable) 


10  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 

PROJECT 

TASK 

ELEMENT  NO 

NO 

NO 

WORK  UNIT 
ACCESSION  NO. 


1 1  TITLE  (Include  Security  Classification) 

DISCRETE  ARMA  MODEL  FOR  NATURAL  RESONANCES  IN  ELECTROMAGNETIC  AND 
ACOUSTIC  SCATTERING 


12  PERSONAL  AUTHOR{S> 

COHEN,  Yuval 


13a  TYPE  OF  REPORT 

Master's  Thesis 


i3b  TIME  COVERED 
FROM _  TO 


15  PAGE  COUNT 

95 


COSAT'  CODES 


GROUP 


SUB-GROUP 


14  DATE  OF  REPORT  {Year,  Month,  Day) 

1990  September 


16  SUPPLEMENTARY  NOTATION  The  views  expressed  in  this  thesis  are  those  of  the 
author  and  do  not  reflect  the  official  policy  or  position  of  the  Depart¬ 
ment  of  Defense  or  the  US  Government. 


18  SUBJECT  TERMS  (Conf/nue  on  reverse  if  necessary  and  identify  by  block  number) 

Natural  Resonances;  Electromagnetic  Scattering; 
Acoustic  Scattering;  Radar  Target  Identifica¬ 
tion;  Prony’s  Method:  ARMA  Model 


19  abstract  {Continue  on  reverse  if  necessary  and  identify  by  block  number) 

Investigations  of  scattered  transient  waveforms  from  conducting  bodies  have 
shown  that  it  is  possible  to  classify  electromagnetic  scatterers.  The  con¬ 
cept  is  based  upon  the  natural  resonance  modes  which  are  part  of  the  scat- 
terer  response  to  an  incident  excitation.  A  new  approach  for  describing 
natural  resonance  modes  using  recursive  systems  is  introduced.  A  discrete 
auto- regressive  moving- average  (ARMA)  type  model  for  the  case  of  the  space- 
time  wave  equation  is  presented.  This  model  results  from  a  finite-differ¬ 
ence  approximation  to  the  wave-equation.  The  ARMA  model  has  spatially- 
independent  coefficients  for  the  temporal  recursive  terms.  Computed 
results  showing  aspect-  and  spatial-independence  of  natural  resonance 
modes,  with  verification  of  the  ARMA  model,  are  also  included.  Applica¬ 
tions  to  target  identification,  using  the  natural  resonant  frequencies  of 
the  target's  echo  signature,  are  considered. 


21  ABSTRACT  SECURITY  CLASSIFICATION 


llJaiNiKKKilMiMil 


20  DISTRIBUTION /AVAILABILITY  0^  ABSTRACT 
IS  UNCLASSIFIED/UNLIMITED  □  SAME  AS  RPT  □  dtiC  USERS 


22a  NAME  OF  RESPONSIBLE  INDIVIDUAL  22b  TELEPHONE  (Include  Area  Code)  I  22c  OFFICE  SYMSOl 

MORGAN,  Michael  A.  408-646-2082 


DD  Form  1473,  JUN  86  Previous  editions  are  obsolete  SECURITY  CLASSIFICATION  OF  this  PAGE 

S/N  0102-LF-014-6603  UNCLASSIFIED 


1 


Approved  for  public  release;  distribution  is  unlimited. 

Discrete  ARMA  Model  for  Natural  Resonances 
in  Electromagnetic  and  Acoustic  Scattering 

by 

Yuval  Cohen 

Lieutenant  Commander,  Israeli  Navy 
B.Sc.  in  Electrical  Engineering,  Technion-Israel 
Institute  of  Technology,  Haifa,  Israel,  1986 

Submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  ELECTRICAL  ENGINEERING 

from  the 


Author: 
Approved  by: 


NAVAL  POSTGRADUATE  SCHOOL 
September  1990 


Yuval  Cohen 

Michael  A.  Morgan,  Th«is  Advisor 

(LAw/yC— _ 

Richard  W.  Adler,  Second  Reader 

_ C,  _ 

Michael  A.  Morgan,  Cnairman 
Department  of  Electrical  and  Computer  Engineering 


ABSTRACT 


Investigations  of  scattered  transient  waveforms  from  conducting  bodies  have 
shown  that  it  is  possible  to  classify  electromagnetic  scatterers.  The  concept  is  based 
upon  the  natural  resonance  modes  which  are  part  of  the  scatterer  response  to  an 
incident  excitation.  A  new  approach  for  describing  natural  resonance  modes  using 
recursive  systems  is  introduced.  A  discrete  auto-regressive  moving-average  (ARMA) 
type  model  for  the  case  of  the  space-time  wave  equation  is  presented.  This  model 
results  from  a  finite-difference  approximation  to  the  wave-equation.  The  ARMA 
model  has  spatially-independent  coefficients  for  the  temporal  recursive  terms. 
Computed  results  showing  aspect-  and  spatial-independence  of  natural  resonance 
modes,  with  verification  of  the  ARMA  model,  are  also  included.  Applications  to 
target  identification,  using  the  natural  resonant  frequencies  of  the  target’s  echo 
signature,  are  considered. 
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1.  INTRODUCTION 


A.  BACKGROUND 

Investigations  of  transient  scattered  waveforms  from  conducting  bodies  have 
shown  that  it  is  possible  to  classify  electromagnetic  scatterers.  Such  research  is 
particularly  applicable  to  inverse  scattering  and  radar  target  identification. 
Conceptual  applications  in  radar  target  identification  have  been  demonstrated  in 
some  laboratories  using  advanced,  high  resolution  radar  techniques.  One  type  of 
technique  is  based  upon  the  natural  resonance  modes  which  are  part  of  the  scattering 
response  to  an  incident  transient  excitation.  In  1971,  Baum  [Ref,  1],  introduced  the 
development  of  the  singularity  expansion  method  (SEM),  which  presents  the  response 
of  a  system  as  a  weighted  expansion  of  complex  natural  modes.  Theoretical  studies 
and  experiments  have  shown  that  these  modes  are  functions  of  the  scatterer  geometry 
and  composition  but  are  independent  of  the  incident  excitation,  including  aspect  and 
polarization. 

The  fact  that  these  natural  modes  are  only  functions  of  the  target  led  to  the 
idea  to  use  them  as  a  data  base  for  the  target  identification  process.  This  concept 
was  introduced  by  Moffatt  and  Mains  in  1974  [Ref.  2].  The  identification  process, 
in  its  most  elementary  form,  includes  a  comparison  against  the  modes  that  have  been 
extracted  and  identified,  using  advanced  signal  processing  techniques,  as  applied  to 
the  target’s  time-domain  scattering  response.  Morgan  showed  that  a  complete 
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description  of  the  scattered  signal,  using  the  conventional  SEM  approach,  is  valid  only 
for  the  "late-time"  portion  of  the  scattered  field  [Ref  3].  This  late-time  scattering 
response  is  due  to  the  source-free  currents  that  remain  after  the  incident  field  has 
completed  its  illumination  of  the  target. 

The  unique  properties  of  the  late-time  response  are  crucial  to  the  development 
of  algorithms  to  identify  the  resonance  modes  of  the  target.  These  natural  resonance 
modes  can  be  represented  by  pairs  of  conjugate  poles  in  the  left-half  complex  s-plane 
of  the  Laplace  transform  domain.  Targets  of  different  geometry',  or  composition, 
have  different  pole  representations.  Several  techniques  have  been  developed  to 
extract  the  dominant  poles  in  the  time  domain  scattering  responses  of  simple  targets. 
Morgan  [Ref  4],  for  example,  introduced  the  classiOcation  of  some  kinds  of 
electromagnetic  scatterers  by  the  annihilation  of  the  natural  modes.  The  advantages 
of  this  technique  over  others  were  achieved,  primarily,  by  using  only  the  late-time 
scattered  field.  The  theory  of  natural  resonance  scattering  is  the  basis  for  this  thesis, 
and  it  is  therefore  explained  in  more  detail  in  Chapter  II. 

B.  PRESENTATION  IN  A  DISCRETE  MODEL 

Based  on  the  theory'  of  the  natural  resonance  scattering  one  can  recognize  the 
late-time  portion  of  the  scattered  signal  as  the  response  of  a  linear  time-invariant 
(LTI)  sy-stem.  These  kinds  of  systems  can  be  numerically  modeled  and  described  by- 
means  of  linear  constant-coefficient  difference  equations  [Ref.  5].  A  discrete  model, 
which  describes  the  source-free  current  distribution,  may  explicitly  present  the  late- 
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time  response  of  the  scatterer  to  an  incident  field  as  a  function  of  its  natural  modes. 
Moreover,  the  currents  may  be  described  by  an  auto-regressive  moving-average 
(ARMA)  type  model,  having  constant  coefficients  for  the  recursive  terms  which 
generate  the  natural  resonance  modes.  This  type  of  model  represents  an  important 
class  of  discrete  systems  which  are  known  as  recursive  systems  since  the  output 
depends  on  previous  values  of  the  output.  Such  a  discrete  model  can  be  employed 
in  tile  process  of  identifying  the  resonance  modes  for  electromagnetic  and  acoustic 
scatterers  of  various  shapes.  The  time-independent  coefficients  of  the  recursive  terms 
in  the  ARMA  model  are  also  the  denominator  polynomial  coefficients  of  the  system 
transfer  function.  Using  the  well-known  concept  of  zeros  and  poles  which  represent 
the  system  frequency  response,  the  roots  of  this  polynomial  are  the  poles  which 
represent  the  complex  natural  frequencies  of  the  resonance  modes  [Ref.  6].  In 
addition,  such  a  discrete  model  gives  a  complete  description  of  the  natural  resonance 
modes  for  a  given  scattering  object  when  limited  by  the  sampling  frequency  imposed 
by  the  Nyquist  theorem.  This  is  the  basic  idea  of  this  thesis,  in  which  such  a  model 
is  set  up  for  the  case  of  the  numerical  finite-difference  solution  of  the  space-time 
wave  equation.  It  is  believed  that  this  is  the  first  time  that  such  an  approach  has 
been  presented. 

C.  THESIS  GOAL 

The  goal  of  this  thesis  was  to  investigate  the  possibilities  of  describing  the 
natural  resonance  modes  of  electromagnetic  and  acoustic  scatterers  by  discrete 
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ARMA  models.  These  ARMA  models  should  have  constant  coefficients  for  the 


recursive  terms  which  determine  the  natural  resonance  modes.  A  model  for  a  one¬ 
dimensional  structure,  such  as  the  damped  string  with  forcing  function,  illustrates  the 
validity  of  the  proposed  approach.  In  addition,  three-dimensional  scattering 
structures  can  be  analyzed  using  a  space-time  finite-difference  method,  which  is  an 
extension  to  the  considered  herein. 

In  this  thesis  an  attempt  have  been  made  to  present  an  overview  of  some 
approaches  for  solving  the  problem  of  the  natural  modes,  and  constructing  the 
required  model.  A  brief  presentation  of  the  theory  of  natural  resonance  modes,  as 
developed  through  the  time-domain  integral  equation,  is  included  in  Chapter  II. 
Chapter  III  describes  the  simplified  numerical  solution,  via  the  vector  potential,  as 
chosen  to  demonstrate  the  possibility  of  using  an  ARMA  model  for  the  case  of  an 
electromagnetic  thin-wire.  This  technique  is  applied  to  the  analogous  problem  of  a 
finite  string.  Chapter  IV  presents  the  ARMA  model  with  several  analytical  and 
computed  results.  Conclusions  and  a  description  of  some  future  questions  are 
discussed  in  Chapter  V. 
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II.  THE  THEORY  OF  NATURAL  RESONANCE  SCATTERING 


A.  INTRODUCTION 

The  general  topic  to  be  considered  in  this  chapter  is  the  natural  resonance 
mode  representation  for  the  induced  current  and  the  scattered  field  transient 
response  of  a  perfectly  conducting  body.  The  case  of  a  thin-wire  is  presented  in  some 
detail.  The  goal  of  this  chapter  is  to  provide  the  theoretical  background  required  for 
constructing  a  discrete  ARMA  model  for  the  natural  resonances  in  electromagnetic 
scattering.  In  addition,  the  features  of  such  a  model  should  agree  with  this  theory. 
A  description  of  time-domain  solutions  relevant  to  this  work  are  included  in  this 
chapter.  There  are  two  main  steps  in  the  process  of  setting  up  the  ARMA  model. 
First,  the  time-domain  numerical  solution  is  required,  and  second,  based  on  the 
numerical  solution,  the  discrete  model  must  be  constructed. 

There  are  two  independent  techniques  available  for  solving  the  problem  of 
transient  scattering.  The  first  involves  the  computation  of  the  frequency-domain 
response  of  the  structure,  followed  by  inverse  Fourier  transformation  to  yield  the 
time-domain  response.  An  alternative  approach  involves  the  direct  formulation  of 
either  partial  differential  or  integral  equations  in  the  time  domain.  One  way  of 
describing  the  scattered  signal  in  terms  of  natural  modes  is  by  using  the  formalism 
known  as  the  Singularity  Expansion  Method  (SEM)  [Ref.  1].  Mittra  and  Van 
Blaricum  showed  that  the  SEM  pole  singularities  of  a  structure  can  be  estimated 


directly  from  its  time-domain  response  [Ref.  7].  In  this  chapter  time-domain  integral 
equations  are  used  general  to  represent  solution  of  the  transient  electromagnetic 
problem.  Also  considered  will  be  the  conceptual  basis  of  mode  representations  for 
fields  and  currents  in  transient  scattering. 


B.  TIME  DOMAIN  INTEGRAL  EQUATIONS 

Consider  the  general  three-dimensional  transient  electromagnetic  scattering 
problem  as  depicted  by  Fig.  1. 


INCIDENT  FIELD 


PERFECTLY  CONDUCTING 
OBJECT 


Figure  1.  Transient  Electromagnetic  Scattering  [After  Ref.  4] 

The  perfectly  conducting  object  is  illuminated  by  an  incident  impulsive  type  plane 
wave  in  free  space.  The  incident  electric  and  magnetic  fields  are  E‘"'  and 
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respectively.  The  integral  equation  approach  represents  the  induced  current  on  the 
surface  of  the  object  in  terms  of  the  incident  field.  There  are  two  fundamental  types 
of  time-domain  integral  equations:  the  Electric  Field  Integral  Equation  (EFIE),  and 
the  Magnetic  Field  Integral  Equation  (MFIE).  Derivations  of  both  EFIE  and  MFIE 
are  described  by  Mittra  in  detail  [Ref.  8]. 

Integral  equation  derivations  begin  with  the  time-dependent  forms  of  Maxwell’s 
equations  in  free  space.  The  continuity  equation  is  used  to  relate  the  current  density, 
Js,  and  the  charge  density  a.  The  scalar  potential,  0,  and  the  vector  potential.  A,  are 
defined  in  terms  of  the  electric  and  magnetic  fields,  respectively.  The  potentials  are 
related  to  each  other  via  the  Lorentz  condition[Ref.  8].  The  wave  equation  can  be 
derived  independently  for  the  scalar  potential  and  for  the  vector  potential  using 
Maxwell’s  equations.  The  sources  of  the  nonhomogeneous  wave  equations  for  the 
vector  and  the  scalar  potentials  are  the  current  density,  Jj,  and  the  charge  density,  a, 
respectively.  The  solution  is  constructed  using  the  Green’s  function,  which  is  the 
impulse  response  in  time  and  space.  The  integral  representation  of  the  time- 
dependent  electric  and  magnetic  fields  are  obtained  by  applying  the  solutions  of  the 
potentials  for  the  general  current  and  charge  distribution. 

The  expressions  of  the  EFIE  and  MFIE,  for  the  scattering  problem  are  finally 
derived  by  applying  the  appropriate  boundary  conditions  of  the  tangential  electric  and 
magnetic  fields  on  the  surface  of  the  perfect  conductor.  In  the  case  of  the  EFIE  the 
sum  of  the  tangential  scattered  electric  field  and  the  incident  electric  field  is  zero. 
The  boundary  condition  in  the  case  of  the  MFIE  is  that  the  total  tangential  magnetic 
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field  is  equal  to  the  surface  current  density.  The  following  equations  are  obtained  for 
the  EFIE  (1)  and  MFIE  (2),  [Ref.  8] 


i^_a_ 

R  dx 


a(r',x)  R 


13  ,  /  X  jR 
-—oir,xy — - 

cR^ 


dx 

O 


ds^ 


(1) 


In 


"■  i  fe 


dx 


(2) 


where  5  is  the  surface  of  the  conductor,  the  unit  vector  which  is  outward  normal  to 
5  is  ft  ,  the  current  density  is  J,,  the  charge  density  is  a,  the  permeability  and 
permittivity  of  free  space  are  /i^  and  respectively.  Further,  the  position  vector 
is  r,  r'  eS  indicates  points  on  the  PEC  surface,  c=  l/(eoMo)^^  js  the  velocity  of  light, 
T=t-R/c  is  the  retarded  time,  and  R=lRl=lr-r'l. 

Both  the  EFIE  and  MFIE  are  Principal  Value  (P.V.)  integrals  because  of  kernel 
singularities  for  r-*r'.  The  special  forms  of  the  time-domain  integral  equations  (1) 
and  (2)  play  fundamental  roles  in  their  solution  construction.  The  main  difference 
between  these  equations  and  their  respective  frequency  domain  equations  is  in  the 
solution  construction.  The  equations  in  frequency  domain  are  handled  numerically 
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by  matrix  inversion,  while  initial-value  techniques  are  applied  to  time  domain  integral 
equations. 

The  solution  of  the  general  integral  equation  for  the  scattering  problem  uses  a 
time-stepping  technique.  Special  cases  such  as  one-  and  two-dimensional  scatterers, 
i.e,  cylindrical  and  wire-type  scatterers,  have  special  forms  of  integral  equations, 
hence  construction  of  their  solutions  differ,  considering  numerical  aspects  and 
accuracy.  For  three-dimensional  scatterers  the  MFIE  is  most  convenient,  while  the 
EFIE  is  used  in  the  case  of  thin  wires  and  thin  surfaces.  This  is  due  to  the  fact  that 
for  solid  surface  structures  the  kernel  in  the  MFIE  is  less  singular  than  the  kernel  in 
the  EFIE.  As  a  consequence,  less  sophisticated  expansion  functions  may  be 
employed  for  representing  the  unknown  current.  On  the  other  nand,  the  MFIE 
becomes  unstable  for  thin  structures.  The  vector  cross  product  between  Jl  and  R  in 
the  Green’s  function  may  lead  to  computational  difficulties  by  virtue  of  the  small 
angle  subtended,  as  in  the  case  of  the  thin-wire  [Ref.  9].  There  is  one  unique  feature 
of  these  integral  equations  for  the  induced  surface  current,  J^,  which  is  important  to 
the  discussion  of  natural  resonance  scattering.  The  unknown  surface  current,  Jj, 
inside  the  integral  in  both  the  MFIE  and  EFIE  has  the  retarded  time  t  argument. 
The  retarded  time  T=t-R/c  is  always  less  then  t  since  the  P.V.  integi.)  'foiuces  the 
point  R=0.  Considering  this  feature,  the  unknown  current  Js(t),  at  any  given  time 
t,  can  be  calculated  from  the  MFIE  in  (2)  using  the  known  incident  field  at  time  t 
combined  with  the  integral  of  terms  which  are  known  from  the  past  histoiy  of  the 
current.  From  another  point  of  view,  the  effect  of  the  source  current  at  any  point  r' 
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is  delayed  by  a  time  R/c  in  affecting  the  current  at  the  observation  point  r.  This 
point  forms  the  basis  of  an  iterative  technique,  termed  dme-stepping,  for  constructing 
the  solution.  The  surface  current  may  be  determined  by  "stepping  on"  in  time,  thus 
eliminating  the  matrix  inversion  required  for  the  numerical  of  the  frequency 

domain  integral  equation.  Figure  2  shows  the  region  o.  ■  i*  'aui  space-time  for 
the  one-dimensional  case  (x,t).  This  region  is  defined  by  ct-l  ^-x’  I  <0  and  denoted 
by  the  shaded  area. 


Ct 


Figure  2.  Region  of  Interaction 

The  same  general  procedure  can  be  similarly  employed  for  solving  the  EFIE 
(1).  The  current  at  the  spatial  point  of  ob.servation,  and  at  some  time  t,  can  be 


in 


calculated  in  terms  of  the  known  incident  electric  field  at  that  same  spatial  point  and 
time,  and  from  the  integral  of  terms  involving  retarded  time  [Ref.  8].  This  point  is 
further  illustrated  in  the  case  of  a  thin-wire. 

C.  NATURAL  RESONANCE  SCATTERING 

The  following  discussion  is  based  upon  the  theory  described  by  Morgan  [Ref. 
4].  The  MFIE  (2)  describes  the  induced  current  on  the  surface  of  the  object  in  terms 
of  the  incident  field  or  "physical  optics"  part,  and  in  terms  of  the  "feedback"  current. 
Figure  3  shows  the  situation  in  transient  scattering.  An  incident  field  with  short  pulse 
duration  illuminates  the  scatterer.  In  the  case  of  a  radar  it  can  usually  be  considered 
a  plane  wave.  Once  illumination  of  the  object  is  completed,  for  t>tL,  H‘"'=0,  and 
only  the  source-free  currents  remain  on  the  object,  generating  the  natural  modes. 
These  modes  are  of  the  form  Jn(r)  exp(s„t),  where  the  natural  resonance  frequencies 
Sn=crj,+jo)n  are  functions  of  the  scatterer  geometry  and  composition.  The  SEM 
representation  of  the  current  takes  the  form  of  a  summation  of  damped  sinusoids 
which  can  be  expressed  using  complex  exponentials  as 

t>h-  P) 

M=— » 

Since  /(r',t)  is  real,  the  complex  exponents  which  are  poles  in  the  complex 
frequency  plane,  and  the  weight  coefficients  ^4^,  come  in  conjugate  pairs.  For  a 
practical  incident  field,  having  a  finite  bandwidth,  only  a  limited  number  (say  N)  of 
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natural  resonance  modes  will  be  "significantly"  excited.  The  scattered  field  which  is 
generated  by  the  induced  current  is  composed  of  two  parts;  an  early-time  driven 
response  and  a  late-time  natural  mode  response  [Ref  4].  The  early-time  scattered 
field  can  be  described  as  a  linear  combination  of  two  terms.  The  first  is  an  aspect- 
dependent  physical  optics  term. 


Figure  3.  SSiort  Pulse  ^yave  liluminntion  [After  Ref.  4] 


The  second  term  describes  the  scattered  field  due  to  the  source-free  current 
distribution.  These  currents  are  integrated  over  the  time-varj'ing  portion  of  the 
surface  to  get  the  contribution  from  all  previously  illuminated  points  which  are  dot- 
shaded  in  Fig.  3.  This  term  can  be  represented  by  the  SEM  expansion  having  the 
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same  exponential  resonance  terms  with  time-varying  coefficients  (SEM  class  2).  The 
late-time  scattered  field  is  due  to  the  source-free  currents  that  remain  after  the 
incident  field  has  completed  its  illumination  of  the  scatterer.  The  late-time  starts 
after  a  delay  of  To=T+2D/c  from  the  initial  impact  of  the  incident  field  on  the 
object,  where  T  is  the  practical  incident  field  pulse  width,  and  D  is  the  length  of  the 
object  in  the  direction  normal  to  the  wavefront  of  the  incident  field.  At  time  T^  the 
integral  is  calculated  over  a  fixed  surface  area,  and  thus  the  SEM  representation  for 
the  late-time  is  a  summation  of  the  same  exponential  resonance  terms  with  constant 
coefficients  (SEM  class  1). 

In  summary,  the  monostatic  transient  scattered  signal  waveform  can  be 
described  in  the  following  form  [Ref.  4] 

>■(0 I  (4) 

where  U(t)  is  the  Heaviside  unit  step  function,  yjj(t)  is  the  early-time  response,  yL(t) 
is  the  late-time  response,  and  N(t)  describes  the  measured  noise  and  other  signal 
pollutants.  The  late-time  portion  of  the  scatterer  response  to  such  an  incident  field, 
with  its  unique  features,  is  considered  in  the  process  of  this  research.  The  discrete 
form  of  the  late-time  response  is  obtained  in  order  to  set  up  the  required  ARMA 
model. 
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D.  THIN  WIRE  CASE 


Consider  the  case  of  a  perfectly  conducting  cylindrical  thin-wire  scatterer.  The 
geometry  of  the  problem  is  shown  in  Fig.  4.  The  radius  of  the  wire,  a,  is  small 
compared  to  the  wire  length,  L.  The  radius,  a,  is  also  small  compared  to  the  shortest 
significant  incident  field  wavelength,  X. 


z 


Figure  4.  Geometry  of  a  Cylindrical  Thin  Wire  Scatterer 

In  tiiis  case  the  thin-wire  approximation  is  applicable;  that  is,  the  azimuthal  surface 
curre-'t  is  negligible  compared  to  the  axially-directed  component.  Then  the  surface 


14 


current  density,  can  be  written  as 


*  27Ifl 


(5) 


where  ^  is  a  tangential  unit  vector  pointing  in  the  axial  direction.  The  current  is  a 
function  of  z  only.  The  incident  electric  field  has  contribution  to  the  surface 
current  only  in  the  2  direction.  Under  these  conditions  the  integral  equation  for  this 
case  can  be  written,  using  the  EFIE  type  (1),  as 


4nE^(z,x)=f-^-^I(z',x)dz' 
i  R  ot 


Uf^f±IXzW)dv'dz' 


^oi  R^  Ldz' " 


1  r  (z-zO  9 


R^  dz'  ' 


lXz!,<)dz! 


(6) 


where  R=r-r^  is  the  vector  from  the  current  element  along  the  z  axis  pointing  to  the 
point  of  observation  on  the  surface  of  the  wire,  lRl=[(z-z’)^+a^]^'^,  while  T=t-R/c 
is  the  retarded  time.  The  integral  in  (6)  is  not  a  P.V.  type  integral  since  x>  *x. 
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There  are  several  approaches  to  solve  the  thin-wire  integral  equation  arising  in 
time-domain  scattering  problems.  The  interpolation  procedure  and  the  finite- 
difference  approach  are  presented  in  [Ref.  8].  The  general  interpolation  procedure 
is  to  subsectionalize  the  thin-wire  by  dividing  it  into  N  segments  and  then  to  define 
a  set  of  basis  functions  for  expressing  the  unknown  current,  I,  in  each  of  these 
segments.  Similar  segmentation  is  also  necessary  for  the  time  domain,  choosing  the 
appropriate  time  interval  with  regard  of  the  spatial  interval.  An  interpolation  scheme 
in  time  and  space  is  then  used  to  express  the  current  at  one  node  (space  and  time) 
in  terms  of  the  current  values  in  the  neighboring  nodes.  The  final  step  is  to  describe 
the  thin-wire  integral  equation  using  the  expansions  and  apply  point  matching  to 
generate  the  desired  matrix  equation. 

A  second  approach  is  to  use  the  finite-difference  method  to  approximate  the 
differential  operators  appearing  in  the  time-domain  integro-differential  equation. 
This  approach  was  introduced  by  Sayre  and  Harrington  [Ref.  10]  where  the  EFIE  is 
written  in  terms  of  the  vector  potential  A.  The  solution  for  the  vector  potential  A  can 
be  based  on  the  finite-difference  method,  as  applied  to  the  driven  wave  equation. 
In  fact,  aside  from  the  specified  boundary  conditions,  the  solution  has  the  same  form 
as  will  be  employed  in  the  acoustic  string  case  in  the  following  chapter. 

£.  SUMMARY 

The  late-time  portion  of  the  scatterer  response  can  be  represented  by  a 
weighted  summation  of  natural  modes.  The  modes  are  functions  of  the  scatterer 
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geometry  and  not  of  the  incident  excitation.  The  EFIE  which  describes  the  surface 
current  on  the  scatterer  can  be  solved  numerically  for  the  case  of  a  thin-wire. 

Although,  in  principle,  the  ARMA  representation  can  be  derived  by  discretizing 
an  integral  equation,  the  full  topological  connectivity  wrought  by  the  Green's 
function  yields  an  imposing  difficulty.  A  better  approach  is  found  by  employing  the 
finite-difference  approach,  with  its  "nearest-neighbor"  connectivity.  This  will  be 
presented  in  Chapter  III. 
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III.  ELECTROMAGNETIC  AND  ACOUS  *  C  SCATTERING  EXAMPLES 

A.  INTRODUCTION 

Open  region  electromagnetic  and  acoustic  scattering  and  radiation  problems  can 
be  formulated  using  the  integral  equation  approach.  Other  techniques  may  also  be 
applied  to  these  problems,  both  in  the  frequency-  and  time-domains.  One  of  these 
methods  is  the  finite-difference  scheme  which  provides  a  convenient  means  for 
deriving  a  time-stepping  algorithm  for  solving  the  EFIE.  In  this  work,  an  exclusive 
use  of  this  method  has  been  utilized  to  formulate*  the  problem  of  transient  scattering 
via  the  time-space  wave  equation,  and  to  demonstrate  the  new  approach  of 
presenting  the  late-time  scattering  response  in  a  discrete  ARMA  model. 

The  objective  of  this  chapter  is  to  describe  the  numerical  solution  chosen  for 
this  research.  The  derivation  of  the  EFIE  expressed  in  terms  of  the  vector  potential 
is  presented.  An  analogous  case  of  a  vibrating  string  of  fixed  length  is  then  described 
in  conjunction  w'th  various  methods  of  solution  to  complete  the  theory  of  the  analog 
form  of  the  problem.  A  discrete  form  of  the  time-space  wave  equation  using  the 
finite-difference  method  is  then  presented  along  with  pertinent  numerical 
considerations. 

The  wave  equation  with  forcing  term  describes  numerous  physical  phenomena 
such  a  driven  finite  string  or  an  illuminated  transmission  line.  A  computer  program 
entitled  TH7.FOR  was  written  to  support  this  research,  and  a  source  listing  is  given 


18 


in  Appendix  A.  The  program  provides  the  solution  of  the  inhomogeneous  time-space 
wave  equation.  The  string  (or  transmission  line)  can  be  excited  by  a  Gaussian  pulse 
from  different  angles  of  incidence.  The  amplitude  and  width  of  the  Gaussian  pulse 
excitation,  as  well  as  the  number  of  segments  on  the  string  can  be  independently 
selected  by  the  user.  Results  in  time-space  are  presented  for  different  cases.  Fast 
Fourier  transforms  (FFT)  are  used  to  obtain  the  frequency-domain  results  from  the 
time- domain  data.  Results  in  these  cases  are  also  presented  showing  agreement  with 
the  basis  of  the  natural  resonance  scattering  theory.  Time-domain  data  provided  by 
the  computer  program  was  also  used  to  check  the  ARMA  model  and  is  presented  in 
Chapter  IV. 

B.  DERIVATION  OF  THE  VECTOR  POTENTIAL  EQUATION 

The  Electr’c  Field  Integral  Equation  (EFIE)  may  be  written  in  terms  of  the 
vector  potential  A.  The  derivation  begins  with  a  slightly  different  form  of  the  EFIE. 
The  scattered  electric  field,  produced  by  the  induced  current,  may  be  written  in  terms 
of  the  scalar  potential  cf),  and  the  vector  potential  A,  as 

£fr.0=~-V4>.  (7) 

at 

The  scalar  potential  (f)  can  be  expressed  in  terms  of  A,  via  the  Lorentz  condition,  as 

V-A(r,f)  +  ® 

at 
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Applying  the  boundary  condition  of  the  electric  field  on  the  surface  of  a  perfectly 
conducting  body  in  order  to  write  Eq.  (7)  in  terms  of  the  scattered  field,  and  using 
the  Lorentz  condition  to  write  Eq.  (7)  in  terms  of  the  vector  potential,  the  tangential 
component  of  Eq.  (7)  can  be  written  as 

t 

=  f rixVV‘Adt'  ,  reS  .  W 

dt 


Equation  (9)  is  differentiated  with  respect  to  time  to  eliminate  the  integral.  This 
yields  the  following  equation 

»  dJE  1  «  A  ^  c  n  m 

nx - =«x - nxVVA,  reS, 

dt  3/2  c2 

where  the  expression  for  the  vector  potential  A  in  terms  of  the  induced  surface 
current  Jj  is 


R 


(11) 


Equation  (10)  may  be  reduced  to  the  one-dimensional  case  of  the  cylindrical  thin  wire 
along  the  x-direction,  with  length  L  and  radius  a,  under  the  assumptions  described 
in  Chapter  II.  The  following  equation  is  obtained 


1  a^A  .  dE^ 
- =  -47tC„ - 

ax^  c"  a/"  "  dt 


(12) 
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The  vector  potential  in  Eq.  (12)  is  defined  in  a  slightly  different  manner  from  Eq. 
(11).  The  vector  potential  A  is  defined  as 


L 


(13) 


with  R=[(x-x')^+a^]^^. 

For  the  purpose  of  this  work  it  was  found  convenient  to  modify  equation  (12). 
The  term  which  includes  the  incident  field  is  replaced  by  the  function  f(x,t)  which 
describes  the  pulse  excitation.  A  term  which  describes  loss  was  added  in  order  to 
investigate  cases  where  damping  was  included  as  well.  The  final  equation,  which  has 
been  solved  numerically  by  computer  program,  has  the  following  form: 


1  , 
dt^ 


(14) 


where  the  function  U(x,t)  satisfies  the  expression,  and  $  is  the  positive  coefficient  of 
the  loss  term.  When  $=0  the  equation  is  reduced  to  the  lossless  case  which  is 
simply  the  nonhomogeneous  wave  equation.  The  last  step  is  to  define  the  boundary 
and  the  initial  conditions  in  order  to  completely  describe  this  problem. 
Homogeneous  boundary  conditions  have  been  set  in  this  case  to  present  a  total 
reflection  at  the  ends  of  a  string.  This  formulation  also  represents  the  cases  of  short 
or  open  circuits  at  the  ends  of  an  excited  TEM  mode  transmission  line.  In  the 
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electromagnetic  case  the  potential  on  the  ends  of  the  single  wire  scatterer  is  non-zero 
and  is  thus  not  properly  represented  by  the  wave  equation  solution  considered  here. 

In  terms  of  the  function  U(x,t)  the  boundary  conditions  may  be  written  as 
U(0,t)=U(L,t)=0.  The  initial  conditions  can  describe  the  voltage  or  current  on  the 
transmission  line  before  the  incident  field  impacts.  There  is  assumed  to  be  no  initial 
voltage  or  current.  Likewise,  the  string  is  initially  at  rest  in  the  acoustic  case.  These 
conditions  are  formulated  by  homogeneous  initial  conditions  as  U(x,0)=0,  and 
9U(x,0)/3t=0,  where  0<  x  <L.  Note  that  homogeneous  initial  conditions  may  be 
replaced  by  some  functions  U(x,0)=h(x)  and  5U(x,0)/9t=g(x),  where  0<  x  <L,  in 
order  to  make  the  partial  differential  equation  (PDE)  homogeneous.  This  point  is 
further  illustrated  by  the  ARMA  model  presentation. 

C.  VIBRATING  STRING  WITH  FORCING  FUNCTION 

The  case  of  a  thin  wire  is  analogous  to  the  problem  of  a  vibrating  string  of  fixed 
length.  This  problem  is  thoroughly  discussed  in  the  literature,  including  [Ref.  10] 
which  covers  this  topic  in  considerable  detail.  The  formulation  of  a  vibrating  string 
problem  is 

±^VM-n  (15) 

with  the  boundary  conditions  U(0,t)=0,  U(L,t)=0,  and  initial  conditions  U(x,0)=h(x), 
and  3U(x,0)/9t=g(x).  Several  methods  are  available  to  solve  the  problem  of  a 


22 


vibrating  string  of  fixed  length.  One  way  to  solve  it  is  by  separation  of  variables. 
This  yields  the  eigenvalues  and  eigenfunctions  corresponding  to  the  homogeneous 
boundary  conditions.  A  Fourier  series  method  is  then  used  to  satisfy  the  initial 
conditions.  This  is  found  to  be  a  convenient  means  to  solve  the  problem  in  many 
cases.  However,  difficulties  may  arise  due  to  the  complexity  of  some  initial  condition 
functions  which  must  be  integrated.  One  applicable  example  in  the  case  of  transient 
scattering  could  be  the  function  which  describes  a  Gaussian  pulse  excitation. 

An  equivalent  technique,  but  in  some  ways  more  useful,  is  the  "method  of 
characteristics".  In  this  approach  the  d’Alembert’s  solution  is  applicable  [Ref.  11], 
with  the  following  form 

x*ct 

U(X.I)  =  X  ±  /  g(x')dx'.  (10 


The  solution  is  valid  for  the  region  defined  by  0<  x-ct  <L  and  0<  x+ct  <L  .  This 
region  is  shown  dot-shaded  in  Fig.  5  which  describes  points  of  position  and  time  such 
that  signals  from  the  boundary  have  not  already  arrived.  The  modification  to  the 
solution  is  made  considering  the  boundary  conditions  which,  in  turn,  imply  multiple 
reflections  as  illustrated  in  Fig.  6.  The  simplest  way  to  obtain  the  solution  is  to 
extend  the  initial  condition  functions  as  odd  functions  (around  x=0)  with  period  2L. 
With  these  odd  periodic  initial  conditions,  the  method  of  characteristics  can  be 
utilized  as  well  as  d’Alembert’s  solution  (16). 
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D.  FINITE-DIFFERENCE  APPROXIMATION  TO  THE  WAVE  EQUATION 


The  discrete  form  of  Eq.  (14)  is  required  for  the  numerical  solution.  In  general, 
numerical  solutions  are  required  in  cases  where  the  shape  of  the  area  of  integration 
or  changes  to  the  boundary  and  initial  conditions  make  analytical  solutions 
impossible.  These  changes  do  not  fundamentally  affect  finite-difference  methods 
although  sometimes  modifications  to  the  methods  are  necessary.  In  this  work  the 
finite-difference  method  was  used  to  approximate  the  derivatives  in  Eq.  (14),  and  to 
make  use  of  the  time-stepping  algorithm  previously  mentioned. 

The  finite-difference  method  may  be  derived  based  on  a  Taylor’s  series 
approach.  The  idea  is  to  approximate  the  function  U(x)  at  a  point  near  x=Xo,  e.g., 
(Xo±  Ax),  using  the  polynomial  approximation  of  the  Taylor  series.  Through  the  use 
of  Taylor  series,  it  is  possible  to  approximate  derivatives  in  various  ways.  A  finite- 
difference  approximation  for  a  derivative  can  be  written  using  forward  difference, 
backward  difference  or  central  difference.  The  central  difference  is  found  to  be  more 
accurate  [Ref.  11]  and  it  was  used  in  this  work. 

Using  the  central  difference  formula,  the  first  partial  derivative  is  written  as 

ax  ^x  ’ 

and  the  second  derivative  as 


,  U(x^+^x,y^)-'^U(x^,yJ+U(x^-Ax,yJ 

- - 


(18) 
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Figure  7  shows  the  discrete  space-time  plane  of  the  finite  string  problem.  The 
string  length  L  is  subdivided  into  N  segments  with  Ax=L/N  being  the  spatial  interval. 
The  time  interval  is  chosen  to  be  At==Ax/c,  where  c  is  the  velocity  of  propagation. 
The  notation  U/,  with  the  spatial  segment  i  and  the  time  segment  j,  is  used  to 
describe  the  function  U(x,t),  where  x=iAx  and  t=jAt. 
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Figure  7.  Space-Time  Discretization 


The  approximated  PDE  for  Eq.  (14)  may  be  written  as 


UU-2UI.UU 


2At 


=//• 


(19) 
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The  relation  between  the  space  interval  and  the  time  interval  may  be  used  in 
Eq.  (19)  to  yield  Eq.  (20).  Equation  (20)  is  written  for  in  a  form  called  the  jrar 
operator.  The  star  operator  is  illustrated  in  Fig.  8. 


Figure  8.  Star  Operator 


Equation  (20)  describes  the  standard  five-point  finite-difference  approximation  to  the 
Laplacian  V"  with  slightly  different  weights. 


where 

A=-2l-,  P=-^-A,  D=-A£ix\  cl=ci6.x  . 
2+cJ  2+cJ 


(20) 


Note  that  when  $ =0,  which  is  the  lo.sslcss  case,  cl=0,  A=l,  P=-l,  and  D=-Ax^.  The 
star  operator  is  used  to  de.scribe  the  finite-difference  equation  in  a  simple  way.  Note 
that  in  this  case  the  value  of  U/  does  not  contribute  since  At=Ax/c.  Given  the  initial 
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and  boundary  conditions,  which  are  homogeneous  in  this  case,  the  values  of  the 
function  U/  at  any  point  (i,j)  can  be  calculated  by  marching  forward  in  time  using  the 
star  operator.  This  is  the  solution  procedure  applied  in  the  computer  program  given 
in  Appendix  A. 

E.  NUMERICAL  CONSIDERATIONS 

In  any  application  of  numerical  solutions,  the  questions  of  accuracy,  stability, 
and  information  bandwidth  must  be  resolved.  This  section  includes  a  discussion  of 
some  numerical  aspects  which  are  relevant  in  this  case. 

The  finite-difference  approximation  to  the  first  derivative  is  consistent,  meaning 
that  the  truncation  error  vanishes  as  Ax-*0.  Hence,  the  exact  solution  is  expected  to 
converge  as  the  number  of  segments  increases.  The  sampling  rates  in  the  spatial  and 
temporal  domains  are  key  factors  in  determining  the  accuracy  of  the  numerical 
solution.  The  spatial  sampling  rate  should  be  high  enough  to  adequately  resolve  the 
spatial  variation  of  the  incident  field  as  it  propagates  past  the  scattered  The  time 
sampling  space  should  be  high  enough  to  adequately  resolve  the  time  variation  of  the 
pulse  excitation.  However,  the  sample  points  in  time  are  not  indeoendent  of  the 
space  interval.  Correlation  between  them  is  required  because  of  equivalence  between 
space  and  time  in  the  retardation  effect.  Since  the  interactions  between  the  currents 
on  different  points  on  the  scatterer  depend  upon  the  velocity  of  propagation,  the  time 
sample  spacing,  At,  must  be  related  to  the  space  sample  interval.  Ax  by  cAt<  Ax  [Ref. 
9].  This  is  also  known  as  the  Courant  stability  condition  for  the  wave  equation  [Ref 
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11].  The  inequality  is  equivalent  to  requiring  that  the  space  sample  points  be  at  least 
as  far  apart  as  the  distance  the  electromagnetic  wave  travels,  with  velocity  c,  in  the 
interval  between  two  sample  points  in  time.  The  relation  between  At  and  Ax 
determine  the  stability  of  the  solution.  The  idea  is  that  the  numerical  process  should 
limit  the  amplification  of  the  initial  conditions.  In  this  simple  case  of  a  finite  string, 
the  time  sample  spacing  is  related  to  the  space  sample  interval  by  cAt=Ax. 

Among  the  factors  which  determine  the  sampling  rate  are  the  shape  and  the 
width  of  the  incident  field  pulse,  the  scatterer  size  relative  to  the  pulse  width,  and  the 
highest  frequency  natural  mode  to  be  obtained  by  the  solution.  A  delta  function 
space-time  impulse  whose  frequency  spectrum  extends  from  zero  to  infinity  with 
uniform  amplitude  is  desired.  However,  this  is  impossible  from  a  practical 
standpoint.  The  approximation  is  made  by  a  Gaussian  impulse  since  it  rapidly  decays 
to  zero.  The  same  property  is  applicable  in  the  frequency-domain,  where  the 
amplitude  rapidly  decreases  with  increasing  frequency.  In  order  to  adequately  resolve 
the  incident  field  in  time  and  space  the  appropriate  sample  spacing,  which  result  in 
reasonably  accurate  numerical  results,  have  been  found  to  be  on  the  order  of  one- 
fifth  and  up  to  one-tenth  the  pulse  width  in  time  and  space,  respectively.  The  pulse 
width  is  determined  by  the  scatterer  size  and  the  highest  frequency  information 
required. 


F.  COMPUTED  RESULTS 


A  computer  program  was  written  to  obtain  time-domain  results  to  support  this 
investigation.  The  program  numerically  solves  Eq.  (14)  using  the  finite-difference 
method.  The  solution  is  determined  using  the  procedure  described  in  previous 
sections  and  the  pulse  shape  is  Gaussian.  Gaussian  pulse  parameters  such  as  pulse 
width,  amplitude,  and  angle  of  incidence  can  be  changed  independently  by  the  user 
to  provide  results  for  different  cases.  Additional  parameters  such  as  string  length, 
number  of  segments  on  the  string,  time  delay  for  peak  pulse  impact,  and  the  number 
of  time  steps  to  be  computed  can  also  be  varied.  The  coefficient  5,  for  the  case  of 
loss,  can  be  altered  by  the  user  within  the  program  to  investigate  lossy  cases.  Note 
that  for  5=0  the  program  will  solve  for  the  lossless  case.  Figure  9  shows  the 
Gaussian  pulse  excitation  for  different  angles  of  incidence.  Broadside  pulse  excitation 
is  shown  in  Fig.  9(a),  and  Fig.  9(b)  shows  the  pulse  excitation  for  a  60  degree  incident 
angle.  In  both  cases  there  are  10  time  samples  in  the  pulse  width,  where  the  pulse 
width  is  defined  between  the  points  at  which  the  amplitude  is  10%  of  the  maximum. 

Figures  10  and  11  show  some  results  in  the  space-time  domain.  Figure  10 
shows  the  displacement  on  the  string  along  200  time  steps  for  the  case  of  broadside 
excitation.  A  one-meter  "electromagnetic  string"  was  assumed,  with  c=3xl0®  m/s. 
The  one-meter  length  was  subdivided  into  15  segments  which  results  in  a  space 
interval  of  1/15  meter.  The  time  interval  is  related  to  the  space  interval  by  the 
velocity  of  jight  ,c,  which  gives  At=(AY/c)=  0.22  ns.  The  pulse  excitation  has  a 
Gaussian  shape  with  pulse  width  of  2.0  ns.  Figure  10a  shows  the  result  for  the 
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(b)  60  Degrees  Excitation 
Figure  9.  Gaussian  Pulse  Excitation 
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lossless  case  and  Fig.  10b  is  the  case  with  loss  added  to  the  system  by  setting  the 
value  of  the  coefficient  $  in  equation  (20)  to  be  non-zero  and  positive.  In  these 
cases,  when  the  string  is  excited  by  a  broadside  Gaussian  pulse,  the  displacement  is 
symmetric  along  the  string.  Symmetry  is  observed  about  the  center  of  the  string  due 
to  the  symmetric  boundary  conditions,  initial  conditions  and  driver.  The  difference 
between  the  lossless  case  and  the  lossy  case  is  in  the  amplitude.  The  displacement 
of  each  segment  is  a  periodic  function  with  a  period  of  2N  time  steps  where  N  is  the 
number  of  segments  on  the  string,  except  in  the  early-time  when  the  incident  pulse 
is  still  exciting  the  string.  This  point  is  illustrated  in  Fig.  11.  In  this  figure,  the 
displacement  of  segments  number  2  and  8  are  presented  for  the  lossy  case  having 
periodic  propagation  except  when  the  string  is  excited  by  the  incident  pulse.  The 
period  is  30  time  steps  which  is  exactly  2N  for  N=15  segments.  The  early-time  is 
about  15  time  steps  including  several  time  steps  before  the  pulse  starts. 

Figures  12a  and  12b  show  results  for  the  case  of  30  degree  incident  angle  on 
the  Gaussian  pulse.  All  other  parameters  are  the  same  as  for  the  case  illustrated  in 
Fig.  10.  As  expected  the  displacement  of  the  string  is  asymmetric.  The  early-time 
in  this  case  is  longer  than  in  the  case  of  broadside  excitation  since  it  takes  more  time 
for  the  incident  pulse  to  complete  its  excitation  of  the  string.  The  time-domain 
characteristics  in  these  results  did  not  change.  Both  have  the  same  period  (2N  time 
steps).  Figure  12a  shows  the  results  for  the  lossless  case  and  Fig.  12b  for  the  lossy 
case. 
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(b)  Lossy  Case 

Figure  10.  String  Displacement  for  Broadside  Excitation 
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(a)  Lossless  Case 


(b)  Lossy  Case 


Figure  12.  String  Displacement  for  30  Degree  Incident  Angle 
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Time-domain  results  were  further  processed  to  obtain  frequency-domain  results 
using  the  Fast  Fourier  transform  (FFT).  The  objective  was  to  check  the  time-domain 
results  against  the  natural  resonance  mode  theory.  The  frequency  spectrum  of  the 
displacement  on  each  segment  was  checked  in  order  to  analyze  the  resonance  mode 
presentation  as  a  function  of  the  position  on  the  string.  The  procedure  includes 
computation  of  the  spectrum  of  the  displacement  as  a  function  of  time  for  each 
segment  on  the  string  using  the  FFT.  A  representative  example  is  shown  in  Figures 
13  and  14,  where  time-domain  data  of  Fig.  12  was  used  in  this  case.  There  are  15 
segments  on  the  string  and  the  incident  angle  is  30  degrees.  Figure  13  shows  the 
frequency  spectrum  of  the  displacement  of  segments  number  3,  6,  and  11.  It  is  clear 
that  the  modes  have  the  same  frequencies  for  these  segments,  although  each  with 
different  energy  depending  on  position.  Figure  14  shows  similar  results  to  that  of  Fig. 
13,  but  for  all  the  segments  on  the  string.  At  this  point  a  tentative  conclusion  may 
be  made  that  the  resonance  mode  frequencies  are  independent  of  the  position  on  the 
string.  Note  that  the  FFT  of  256  points  was  taken  using  the  data  including  the  short 
early-time  data.  In  this  case,  the  early-time  data  may  be  used  since  this  portion  also 
includes  information  about  the  modes.  Figure  15  shows  the  frequency  spectrum  in 
the  case  of  broadside  excitation  and  with  the  same  parameters  as  in  the  case  of  30 
degree  incident  angle  (Fig.  13).  The  differences  between  the  results  are  explicitly 
presented.  In  the  case  of  broadside  excitation,  only  even  modes  are  excited  while  in 
the  case  of  the  30  degree  incident  angle  both  even  and  odd  modes  are  excited. 
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Relative  Power 


Figure  13.  Resonant  Frequencies  (30  Degrees) 
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Relative  Power 


Figure  15.  Modes  or  All  Segments  (Broadside) 
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The  fact  that  these  modes  are  aspect-independent  may  be  checked  in  the 
frequency-domain  by  using  the  time-domain  data  for  different  incident  angles. 
Figures  16  and  17  show  results  for  various  incident  angles.  The  string  was  subdivided 
into  N=25  segments.  The  Gaussian  pulse  width  was  set  to  2  nanoseconds  in  order 
to  allow  15  samples  in  the  10%  pulse  width.  These  parameters  result  in  a  sampling 
frequency  to  be  7.5  Ghz,  thus  providing  a  frequency  limit  of  3.75  Ghz,  as  determined 
by  the  Nyquist  sampling  theorem.  However,  the  incident  pulse  bandwidth  is  less  than 
1  Ghz,  hence  the  modes  are  presented  only  within  this  bandwidth.  The  string  was 
excited  from  incident  angles  of  0,  25,  50,  and  75  degrees.  Figures  16  and  17  show 
results  of  a  256  point  FFT  of  the  time-domain  data  for  the  same  segment  in  each 
figure.  These  figures  show  results  obtained  on  segment  number  22  and  6, 
respectively.  In  both  cases  the  modes  excited  for  different  incident  angles  appear  to 
be  the  same.  Note  that  for  the  0  degree  incident  angle  (solid  line),  only  the  even 
modes  are  excited. 
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Figure  16.  Modes  on  Segment  22 


4J 


Frequency  (hz) 


Figure  17.  Modes  on  Segment  6 


rv.  THE  DISCRETE  ARMA  MODEL 


A.  INTRODUCTION 

Based  on  the  theory  of  natural  resonance  scattering,  a  radar  target  can  be 
considered  as  a  Linear  Time-Invariant  (LTI)  system.  There  are  several  ways  of 
describing  LTI  systems.  A  linear  system  can  be  described  by  its  impulse  response, 
by  means  of  linear  difference  equations,  a  system  diagram  or  by  the  system  transfer 
function  with  poles  and  zeros.  The  system  transfer  function  model  is  extensively  used 
to  present  natural  resonance  modes  in  transient  scattering  research.  As  previously 
mentioned,  these  modes  are  functions  of  the  scatterer  geometry  and  composition,  and 
are  independent  of  the  incident  excitation.  A  finite  subset  of  the  scatterer  poles  can 
possibly  be  used  to  represent  the  scatterer  in  the  process  of  discrimination. 

An  alternative  linear  system  model  may  be  used  to  describe  a  scatterer  in  order 
to  find  its  natural  modes  in  a  relatively  simple  way.  This  new  approach  describes  a 
system  by  means  of  difference  equations.  Different  scatterers,  which  are  considered 
LTI  systems,  may  be  described  by  means  of  linear  constant-coefficient  difference 
equations.  Once  such  a  model  has  been  set  up  for  a  specific  scatterer,  the  natural 
modes  can  be  directly  determined  by  the  coefficients  of  the  differential  equation. 

In  this  work  such  a  model  has  been  set  up  for  the  late-time  response  of  a 
vibrating  string  with  forcing  function.  An  alternate  physical  problem,  which  fits  the 
mathematical  model,  is  an  illuminated  TEM  mode  transmission  line  with  either  open 
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or  shorted  ends.  In  the  equivalent  electromagnetic  case  the  late-time  starts  after  the 
incident  field  has  completed  its  illumination  of  the  target,  and  is  completely  described 
using  a  weighted  expansion  of  complex  natural  modes.  During  this  period  the 
scatterer  acts  as  a  recursive  system,  generating  the  natural  resonance  modes.  The 
output  of  this  system  depends  only  on  previous  values  of  the  output.  This  type  of 
system  can  be  described  by  an  auto-regressive  moving-average  (ARMA)  model.  The 
coefficients  for  the  recursive  terms  are  constant  since  the  feedback  mechanisris  on 
the  scatterer  are  constant  and  time-independent.  Moreover,  the  ARMA  model  which 
describes  the  system  has  the  same  recursive  coefficients  for  all  spatial  points,  as 
expected. 

Two  cases  are  considered:  lossless  and  lossy  .  The  analytical  solution  and  the 
numerical  solution  for  the  natural  modes  in  both  cases  are  presented  including 
comparisons  with  computed  results.  Development  of  the  ARMA  model  for  each  case 
is  then  presented.  Demonstrations  of  these  developments  are  included  through 
examples  for  both  the  lossless  and  lossy  cases. 

A  variety  of  methods  are  currently  used  to  estimate  ARMA  parameters.  These 
methods  have  been  used  to  estimates  poles  of  scatterers  by  applying  them  to  given 
data  obtained  by  measuring  the  backscattering  signal.  In  this  work  only  the  basic 
version  of  the  Prony’s  method  is  applied  to  data  obtained  by  the  numerical  solution 
described  in  Chapter  III.  The  objective  was  to  verify  the  development  of  the  ARMA 
model.  Results  of  the  ARMA  model  are  compared  with  results  obtained  using  the 
Prony’s  method  for  pole  estimation  [Ref.  13). 
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B.  THE  PROBLEM 


The  problems  of  transient  electromagnetic  scattering  or  the  equivalent  problem 
of  acoustic  scattering  may  be  treated  in  the  sampled  signal  case  by  considering  the 
scatterer  as  a  linear  system  whose  input  and  output  satisfy  a  linear  constant- 
coefficient  difference  equation  of  the  form 

*=1  i.O 

For  transient  electromagnetic  scattering,  it  is  assumed  that  no  surface  current  is  on 
the  scatterer  before  the  scatterer  is  excited  by  the  incident  field,  while  the  scatterer 
is  initially  at  rest  in  the  acoustic  case.  Therefore,  the  system  may  be  considered 
casual,  linear,  and  time-invariant.  The  a's  and  the  b's  in  this  case  are  real  constants 
and  the  difference  equation  in  Eq.  (21)  can  be  used  to  compute  the  output 
recursively  [Ref.  6].  Considering  the  late-time  case,  the  input  x(n)  is  zero  for  n>no, 
where  the  discrete  time  n^,  corresponds  to  T^  in  the  analog  case.  This  is  the  time  at 
'vhich  the  incident  field  has  completed  its  illumination  of  the  scatterer.  Hence,  the 
difference  equation  for  the  late-time  has  the  form 

N 

y(.nhY^a^y(,n-k)  ,  n>n^.  (22) 

*=i 

Equation  (22)  describes  the  system  model  for  the  late-time  where  the  unknowns  are 
the  coefficients  Aj,  a^,  and  the  number  of  delays  N.  Figure  18  shows  the  form 
of  the  system  diagram,  which  serves  as  a  graphical  way  of  representing  the  same 
information  contained  in  the  difference  equation  (22).  Equation  (22)  is  referred  to 
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as  the  homogeneous  part  of  the  general  form  of  an  Nth-order  linear  constant- 
coefficient  difference  described  by  Eq.  (21)  with  a^=\. 


Figure  18.  Realization  of  the  General  Model  [After  Ref.  5] 

The  homogeneous  equation  (22)  has  a  family  of  solutions  of  the  form 

)>)«i 

where  are  complex  numbers.  A  unique  solution  requires  a  set  of  N  auxiliary 
conditions  since  th°re  are  N  undetermined  coefficients.  Substituting  Eq.  (23)  into  Eq. 
(22),  the  complex  lumbers  z^,  must  be  roots  of  the  polynomial 

Eo.r'-o,  (24) 

*=0 

assuming  that  all  N  roots  of  the  polynomial  in  Eq.  (24)  are  distinct.  Based  on  the 
theory  of  natural  resonance  scattering,  these  roots  are  the  poles  in  the  z-transform 
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as  tlie  homogeneous  part  of  the  general  form  of  an  Nth-order  linear  constant- 
coefficient  difference  described  by  Eq.  (21)  with  0^=1. 


Figure  18.  Realization  of  the  General  Model  [After  Ref.  5] 

The  homogeneous  equation  (22)  has  a  family  of  solutions  of  the  form 

W«1 

where  are  complex  numbers.  A  unique  solution  requires  a  set  of  N  auxiliary 
conditions  since  tH°re  are  N  undetermined  coefficients.  Substituting  Eq.  (23)  into  Eq. 
(22),  tlie  complex  lumbers  z^,  must  be  roots  of  the  polynomial 

Ea*z-*=0,  (24) 

*=0 

assuming  that  all  N  roots  of  the  polynomial  in  Eq.  (24)  are  distinct.  Based  on  the 
theory  of  natural  resonance  scattering,  these  roots  are  the  poles  in  the  z-transform 
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X  and  t  may  be  separated  by  substituting  the  solution  (26)  into  the  PDE  (25)  to 
obtain  the  eigenvalue  problem 


dx^ 


(27) 


with  the  boundary  conditions  Wn(0)=«„(L)=0.  The  general  solution  for  the 
eigenvalue  problem  has  the  form 

u„(x)=Cye^'^‘+C2e'^'^‘. 

The  boundary  conditions  are  applied  to  obtain  the  solutions  for  s^  and  the  associated 
modes.  For  the  lossless  case  Sn=j(i)„,  where  Wn=n7rc/L  for  n=±l,  ±2,  ....  Natural 
resonance  modes  have  the  form 


U„ix,t)=sm(!fx)e’'^^\  (29) 

The  trivial  solution  is  obtained  for  n=0.  The  solution  for  the  initial  value  problem 
may  be  obtained  by  writing  the  final  solution  as  a  superposition  of  natural  modes 

«e 

E  .  n*0.  (^C) 

The  constants  ^4^  which  are  usually  complex,  can  be  found  by  applying  the  initial 
conditions  and  representing  them  by  the  appropriate  Fourier  series,  then  matching 
term  by  term. 
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An  alternative  way  to  find  the  natural  mode  expansion  for  the  initial  value 
problem  is  by  using  the  fact  that  these  modes  appear  in  conjugate  pairs  for  real  initial 


values.  In  this  case  the  solution  can  be  represented  by 

m 

U(x,t)  =  Yi  c„sin(^x)cos(o>^/+0„), 

/i=i  X. 

where  U.^=U^.  The  real  coefficients  Cf^  and  the  unknown  phase  (f)„  are  determined 
by  the  initial  conditions  using  Fourier  series  to  represent  the  initial  condition 
functions  and  matching  term  by  term. 

For  the  lossy  case  the  solution  represented  by  Eq.  (26)  may  be  applied, 
Mn(x)=sin(n7rx/L),  but  a  slightly  different  solution  is  obtained  for  the  modes.  This 
solution  is  substituted  into  the  lossy  wave  equation  which  has  the  form 

^.±iK-lK-o  (32) 

dx^  dt  ’ 


where  $>0  is  the  coefficient  for  he  loss  term.  The  eigenvalue  problem  for  the 
spatial  variable  x  has  a  different  form,  yielding  a  complex  set  of  solutions  for  s^. 
Substituting  Eq.  (26)  into  Eq.  (32)  to  separate  the  variables  x  and  t,  the  following 
equation  is  obtained  for  the  variable  x: 


^ujx) 

dx^ 


(33) 
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The  solution  for  s„  is  found  as  roots  of  the  quadratic  equation 

(34) 

D.  NUMERICAL  SOLUTION  FOR  RESONANCE  MODES:  LOSSLESS  CASE 
The  solution  for  the  natural  resonance  modes  involves  the  finite  difference 
discrete  form  of  the  wave  equation,  using  the  star  operator.  Figure  19  shows  the 
space-time  discrete  domain  as  applied  to  this  development.  The  finite  difference 
approximation  results  in  a  discrete  equation  of  evolution  which  may  be  written  for  the 
lossless  case  in  the  following  form 

U(k)  =A‘U(k- 1)  -  U(k-2) ,  (35) 

where 


UiX^ytf) 

0  1  0  •••  o' 

10  10  0 

II 

f 

II 

0  •••  •••  0 

0  0  10  1 

p  0  1  p 

The  vector  U(k)  is  composed  of  unknown  nodal  values  at  the  k-th  time  step  and  A 
is  an  M  X  M  sparse  matrix  with  ones  along  the  two  diagonals  adjacent  to  the  primary 
diagonal,  regardless  of  size.  The  number  of  unknown  nodal  values  is  M  which  is 
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related  to  the  number  of  segments  by  M=N'-1,  where  N'  is  the  number  of 
segments.  The  discrete  form  of  the  natural  mode  solution  has  the  separable  form 


ft 


where  C/„  is  a  constant  in  k,  and  s„=j<On  in  the  lossless  case. 


(36) 


I- 


Figure  19.  Space-Time  Discretization  for  Mode  Solution 
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Both  w„  and  should  have  acceptable  values  to  solve  Eq.  (35).  The  solution  for  the 
natural  inodes  may  be  obtained  by  substituting  the  solution  in  Eq.  (36)  into  Eq.  (35) 
to  yield  the  following  eigenvalue  problem 


(37) 

Equation  (37)  represents  M  linearly-independent  equations  for  the  M  unknown 
eigenvalues  A-n.  After  solving  Eq.  (37)  for  the  eigenvalues,  the  modes  can  be  solved 
by  the  relation  between  and  A.^,  via  Eq.  (37),  as 


2 

n  It  n 


(38) 


Equation  (38)  solves  the  natural  modes.  The  number  of  modes  is  2M  since  for  every 
value  of  A-n  there  are  two  solutions  for  the  Z„  which  are  conjugate  values.  This  can 
be  explained  by  the  fact  that  in  the  finite  difference  equation  (35)  there  are  2M 
degrees  of  freedom  since  there  are  M  unknowns,  and  each  unknown  requires  two 
previous  values.  From  this,  it  is  concluded  that  2M  coefficients  are  expected  to 
appear  in  the  ARMA  model.  In  terms  of  number  of  segments,  the  number  of 
coefficients  in  the  ARMA  model  is  given  by  N=2(N'-1),  where  N'  is  the  number  of 
segments. 

The  initial  condition  solution  may  be  obtained  by  a  superposition  of  the  modes 
as 


(39) 
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where  N  is  the  number  of  modes.  The  A„  terms  are  excitation  dependent  amplitudes 
and  Mn(k)  are  excitation  independent  modes.  There  must  be  at  least  N=2M  modes 
to  allow  the  superposition  (39)  to  be  complete  for  any  given  initial  conditions, 

E.  NUMERICAL  SOLUTION  FOR  RESONANCE  MODES:  LOSSY  CASE 
The  same  solution  procedure  followed  in  the  lossless  case  is  applied  to  the  lossy 
case.  The  finite  difference  equation  may  be  written  as 

TO  =  Cj4*TO- 1)  -  C2TO-2) ,  (40) 


where  Ci=2/(2+A,c5Ax),  C2=-Cj+(c5Ax)/(2+c5Ax),  and  both  A  and  U{k)  have  the 
same  form  as  for  the  lossless  case.  Note  that  for  the  lossless  case  (5=0)  Ci=l  and 
C2=-L  The  same  solution  (36)  is  applied  with  Sn=crn +)<!)„.  By  substituting  the 
solution  (36)  into  the  PDE  (40)  and  separating  the  variables,  the  following  eigenvalue 
problem  is  obtained 


Equation  (41)  represents  M  linearly-independent  equations  for  M  unknown 
eigenvalues.  Solving  Eq.  (41)  for  the  eigenvalues,  the  modes  can  be  obtained  by 
using  the  relation  between  and  in  Eq.  (41)  as 


(42) 
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The  initial  value  problem  may  be  solved  by  a  superposition  of  the  natural 
modes  as  in  the  lossless  case,  with  N=2M  modes  to  allow  the  superposition  (39)  to 
be  completed  for  any  given  initial  conditions. 


F.  ARMA  MODEL  DEVELOPMENT 

The  numerical  solution  shows  that  the  modes  are  spatially-independent,  i.e,  all 
nodes  have  the  same  Z^.  Hence  the  ARMA  model  may  have  the  following  form 


N=:2M 


m-l 


(43) 


where  the  a^'s  are  the  same  for  all  spatial  nodes,  even  those  next  to  the  boundaries, 
and  U(k)  is  the  superposition  of  all  Un-Zn*',  for  ne[-M,M].  Note  that  the  number  of 
the  coefficients  is  N=2M,  where  M  is  the  number  of  spatial  nodes  along  the  string 
excluding  those  on  the  boundaries.  The  unknown  coefficients  can  be  determined  by 
applying  the  z-transform  to  the  recursion  equation  (43),  which  yields 

msl 


Equation  (44)  may  be  written  for  the  coefficients  in  the  form 


N 


ZN  v-v  ,yN-m 


(45) 


ffici 
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or  in  the  form 


^  ^N-\  ,  ^N-i  „  r  „  _n 


(46) 


Factoring  the  left  hand  side  of  Eq.  (46)  into  first-order  terms  gives 


(Z-Zi)(Z-Z.i)(Z-Z2)(Z-Z_2)-(Z-Zj^)(Z-Z.jj,) 


(47) 


where  Z„  are  the  poles  and  Z„=Z.„*.  The  last  step  is  to  use  the  results  of  the  natural 
modes  in  Eq.  (47)  generating  the  polynomial  and  comparing  term  by  term  to  Eq.  (46) 
to  find  the  values  of  the  coefficients. 

The  validation  of  the  ARMA  model  can  be  obtained  by  showing  that  Eq.  (43) 
can  be  derived  directly  from  the  equation  of  evolution  (35)  by  working  backwards. 
This  point  is  further  demonstrated  using  an  example  for  the  lossless  case. 


G.  VALIDATION  EXAMPLE:  LOSSLESS  CASE 

Consider  the  case  of  N'  =  4  segments  of  an  undamped  string.  The  number  of 
unknowns  is  M=N'-1=3.  The  vector  of  the  unknown  nodal  values  at  the  k-th  time 
step  has  the  form 


m= 


H(x,,fj) 

U(X2,t^) 

UiXytD 
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The  equation  of  evolution  in  this  case  has  the  form 


with  the  matrix 


U{k)=A'U{k-\)-UQc-2), 


0  1  0 
1  0  1 
0  1  0 


(48) 


The  number  of  discrete  natural  modes  is  N=2M=6.  The  eigenvalue  problem  with 
the  form  of  Eq.  (37)  is  obtained  by  substituting  the  separable  form  (36)  into  the 
equation  of  evolution  (48).  The  eigenvalue  problem  has  the  following  form  after 
rearranging  the  equation 


(49) 


Substituting  A  and  solving  for  the  eigenvalues  gives  three  solutions:  A.n=0,  and 
X^=±'/2.  The  solution  for  the  modes  is  obtained  via 

MO.  A  -4- 

The  solution  has  six  modes:  e*^^^'’.  Figure  20  shows  the  six  poles 

in  the  z-plane.  The  poles  are  on  the  unit  circle  since  this  system  is  lossless.  Since 
Zn=ej®,  where  6=&)„At,  the  modes  are  the  same  as  in  the  exact  case 
«„ = n7r/4  A  t = nTTC/L. 
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The  ARMA  model  has  the  form 


6  _ 


m^l 


The  z-transform  of  Eq.  (51)  is 


(51) 


(v-«) 


msi 


The  coefficients  a„,  can  be  found  by  factoring 


(52) 


(Z-Z^KZ-ZlKZ-Z^XZ-Z^KZ-Z^iZ-Z;),  (53) 

with  the  known  solutions  for  the  Z^.  The  polynomial  obtained  by  this  procedure  is 

Z®+Z^+Z^+l,  (5^*) 

which  gives  the  coefficients:  fl,=0,  02=-!,  0^=0,  04=-!,  ^5=0,  and  a^=-l.  The 
ARMA  model  for  this  case  has  the  form 


U{k)  =  -  Uik-2)  -  Uik-4)  -  U{k-6) .  (55) 

Verification  of  the  ARMA  model  may  be  done  by  using  the  finite-difference 
equation  of  evolution  (48)  to  find  the  coefficients  in  Eq.  (55).  The  procedure  is  to 
use  Eq.  (48)  to  write  f/(k-l),  U(k-3),  and  U{k-5)  to  finally  obtain  the  form  in  (55). 
This  is  shown  in  Appendix  B. 

In  a  case  of  a  5-segment  string,  there  are  4  unknowns  for  each  time  step,  and 
therefore  8  modes.  The  poles  obtained  in  this  case  are  Z„=e.~^~'’^,  e*^^,  and 

g*j4r/5  poles  are  shown  in  Fig.  (21)  on  the  unit  circle  of  the  z-plane. 
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H.  VALIDATION  EXAMPLE:  LOSSY  CASE 


Consider  the  case  of  N'=4  segments  on  a  damped  string.  The  number  of 
unknowns  is  M=N^-1=3.  The  vector  of  the  unknown  nodal  values  at  the  k-th  time 
step  has  the  form  of  that  in  the  example  of  the  lossless  case.  The  equation  of 
evolution  in  this  case  has  the  form 


where 


UQc)  =  c^A-Uik- 1)  -  Cj  U{k-1) , 


2 

Ci  = - ,  c,  =  — - c, . 

2+clAx  1+c^Ax 


(56) 


The  same  matrix  .4  as  in  Eq.  (48)  is  applicable  in  this  case.  The  number  of  discrete 
natural  modes  is  N=2M=6.  The  eigenvalue  problem  with  the  form  (37)  is  obtained 
by  substituting  the  separable  form  (36)  into  the  equation  of  evolution  (56).  The 
eigenvalue  problem  has  the  same  form  as  in  the  lossless  case.  Substituting  A  and 
solving  for  the  eigenvalues  gives  three  solutions:  A.n=0,  and  k„=±c^V2.  The  solution 
for  the  modes  is  obtained  via 


~^n  Cyy/2. 

The  solution  has  six  poles, 


(57) 


, _  *7-7  ±\/2c,  ±i/2c?+4c. 
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Figure  22  shows  these  six  poles  in  the  z-plane.  The  poles  are  inside  the  unit  circle 
and  on  a  circle  representing  equal  loss  for  all  modes.  In  the  electromagnetic  case, 
as  the  frequency  increases  the  higher  frequency  modes  have  more  loss. 


Z-Plane 


Real  z 

Figure  22.  Modes  of  4  Segment  Damped  String 
The  ARMA  model  has  the  form 

6 

m>l 

The  same  procedure  used  to  find  the  coefficients  for  the  lossless  case  is  applied  in 
the  lossy  case.  The  coefficients  can  be  found  by  factoring  the  solutions  for  Z^. 
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The  coefficients  in  this  case  are:  ai=0,  fl2=3c2+2cj^,  ^3=0,  fl4=-3c2^+2c2Ci^,  ^5=0, 
and  Note  that  the  non-zero  o’s  are  negative  numbers  with  decreasing 

absolute  values  as  their  index  increases.  The  ARMA  model  in  this  case  has  the  form 

U{k)  =  (3c2+2cf)  U(k-1)  -Ocl+lc^c^  Uik-4)  +  clU(k-6) .  (^9) 

Checking  the  coefficients  by  using  Ci=l,  and  C2=-l  yields  the  same  coefficients 
obtained  in  the  lossless  case. 

1.  COMPUTED  RESULTS 

The  objective  of  this  part  of  the  work  was  to  obtain  additional  verification  for 
the  ARMA  model.  Ttie  idea  was  to  apply  an  algorithm  which  can  estimate  ARMA 
model  coefficients  to  the  time  domain  data  generated  by  the  computer  program 
TH7.FOR.  There  are  several  methods  available  to  estimate  these  parameters. 
Among  them  is  the  Prony’s  method  [Ref.  13].  The  Prony’s  method  implemented  via 
a  computer  program  entitled  TEST.M.  A  source  listing  is  given  in  Appendix  C.  The 
program  implements  the  following  procedure  for  each  segment  on  the  wire. 

Using  the  time-domain  results  of  a  given  segment,  the  following  M  x  M  matrix 
and  vector  are  generated 


KCt-l) 

u(t-2)  . 

u(t-M) 

u(t) 

?= 

u(t-2) 

uit-3)  ■ 

,  v= 

uit-l) 

u(,i-M-l)  ■ 

-  u(t-2M-l) 
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where  M  is  the  number  of  coefficients,  and  t  is  the  time  step  at  which  the  ARMA 
model  is  applied.  The  time-step  t  should  be  "late"  enough  so  the  time  step  t-2M-l 
is  in  the  late-time  portion  of  the  time-domain  solution,  u(t)  of  the  given  segment.  The 
vector  a  is  composed  of  the  coefficients  Oj,  fl2>  The  coefficients  are 

computed  by 

a=p-^'V, 

where  P‘^  is  the  inverse  matrix  of  P.  This  method  was  applied  to  various  cases.  The 
results  show  that  the  coefficients  are  the  same  for  all  segments.  The  coefficients  are 
also  the  same  for  any  time-step  t  in  the  late-time,  t-2M-l>no,  where  n^  is  the  initial 
time-step  of  the  late-time  portion  of  the  system  response.  For  all  lossless  cases  the 
results  show  coefficients  with  the  form  ai=0,  3,=-!,  33=0,  84=-!,  ...a^j=-l.  In  all 
lossy  cases  the  results  are  a,=0  for  odd  i,  and  decreasing  absolute  values  of  negative 
numbers  a„  for  even  i.  In  some  cases  the  results  were  different  from  segment  to 
segment.  In  those  cases  where  the  results  produced  segment-dependent  coefficients 
the  frequency  response  showed  that  not  all  of  the  modes  where  excited.  A 
representative  result  is  shown  in  Fig.  23.  In  this  case,  the  damped  string  was 
subdivided  into  11  segments,  hence  the  number  of  modes  is  20.  The  same 
coefficients  were  obtained  for  all  segments,  for  various  incident  angles,  and  at  each 
time-step  in  the  late-time. 
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VI.  CONCLUSIONS 


The  solution  of  electromagnetic  natural  resonance  modes  by  means  of  finite 
difference  equations  promises  to  be  useful  and  effective.  The  ability  to  construct 
discrete  ARMA  type  models  of  complex  scatterers  may  be  applicable  to  estimating 
their  electromagnetic  signature.  In  this  thesis  a  simple  ARMA  model  for  the  case  of 
the  one-dimensional  wave  equation  was  investigated  in  order  to  examine  this  new 
approach.  Complex  three-dimensional  scattering  structures  can  be  treated  using  a 
space-time  finite-difference  approach  which  is  an  extension  to  that  investigated 
herein.  The  ARMA  model  was  constructed  for  the  late-time  response  when  the 
scatterer  acts  as  a  recursive  LTI  system.  The  process  of  constructing  such  a  model 
included  iv/o  main  steps.  First  the  scatterer  should  be  presented  in  a  discrete  form. 
In  this  work  the  finite-  difference  approach  was  presented  as  well  as  the  time-domain 
integral  equation.  The  finite-difference  approach  is  found  to  be  convenient  for  the 
construction  of  the  required  model.  The  advantage  of  this  method  is  in  the  locally- 
connected  discrete  form  which  can  be  effectively  used  in  the  next  step.  In  the  second 
step,  an  algorithm  is  constructed  to  explicitly  present  the  natural  modes.  In  this 
work,  the  algorithm  is  based  upon  the  star  operator  obtained  from  the  finite 
difference  approach.  It  was  shown  that  the  star  operator  was  used  recursively  to 
yield  the  ARMA  model,  with  constant  coefficients. 
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Results  of  this  work  have  shown  again,  in  the  case  of  acoustic  scattering,  that 
natural  resonance  modes  of  a  scatterer  are  spatial-  and  aspect-independent.  The 
scatterer  acts  as  a  complete  recursive  system  only  in  the  late-time.  Since  the  early- 
time  is  important  because  of  practical  considerations  of  signal  to  noise  ratio,  further 
work  is  required  to  extend  the  model  to  and  to  examine  its  complexity  in  the  early- 
time.  The  computer  program  entitled  TH7.FOR  can  be  modified  to  help  this 
investigation. 

The  electromagnetic  thin-wire  scattering  case  was  formulated  via  the  vector 
potential.  However,  the  potential  on  the  ends  of  the  wire  is  unknown.  Assuming 
homogeneous  boundaiy  conditions  reduces  the  problem  to  that  of  a  vibrating  finite 
string  with  zero  displacement  at  the  ends.  Extension  to  the  electromagnetic  case  may 
be  made  by  attempting  to  discretize  the  current  on  the  wire  using  the  EFIE.  In  the 
case  of  currents,  the  boundary'  conditions  are  zero. 

A  better  approach  is  found  by  employing  the  finite-difference  approach.  In  this 
method  the  spatial  point  current  is  given  in  terms  of  its  "nearest-neighbor"  currents 
yielding  the  applicable  star  operator.  The  ARMA  model  is  constructed  using  the  star 
operator. 
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APPENDIX  A.  SPACE-TIME  WAVE  EQUATION  PROGRAM 


The  program  entitled  TH7.FOR  numerically  solves  the  displacement  of  a  finite 
string.  The  equation  for  the  position  f/(x,t)  obeys  the  wave  equation  including  a  term 
of  loss 

dx^  dt^  dt  ^ 

The  program  uses  the  finite-difference  method,  with  central  difference,  to  solve 
the  partial  differential  equation  (PDE).  The  string  with  length  L  meters  is  subdivided 
into  N  segments.  The  boundar'  conditions  f/(x=0,t)  and  l7(x=L,t)  are  set  to  zero. 
The  excitation  is  a  Gaussian  pulse  with  the  form 

The  pulse  width  is  set  between  "10%"  points  (T1  in  the  program),  where  Theta  is  the 
incident  angle.  The  time-delay  to  peak  pulse  impact  is  T2.  The  times  Tl,  T2  and 
Theta  can  be  defined  in  the  program,  along  with  the  amplitude  of  the  Gaussian  pulse, 
GA.  Both  space  and  time  intervals  are  defined  as  S  and  T,  respectively. 

c  VARIABLES  DECLARATION 
c 

INTEGER  N.NNEXT.Th 

REAL  L,S,GA.C1,R,A.B,P.D,U1.U2,Tmax,T1,T2,A2.A1,PL(403) 

REAL  PI,Tht,T3 
CHARACTER*64  TITLE 
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DIMENSION  U(1 03,403),  F(1 03,403) 
c 

write(*,*) 
write(*,*)’ 
write(*,*) 
write(*,*)’ 
write(*,*) 
v/rite(*,*)’ 
write(*,*)' 
write(*,*)’ 
write(*,*)’ 
write(*,*)’ 
write(*,*)’ 
write(*,*V 
write(*,*)’ 
write{*,*)’ 
write(*,*)’ 

Wfite(*,*)’ 
write(*,*)’ 
wrrte(*,*)’ 
write{*,*) 
write(*,*)' 
write(*,*) 
write(*,*) 
pause 
c 

c  INITIALIZATION 
c 

WRITE{*,*) 

WRITE(*,*)  ’Enter  number  of  your  MONITOR  Type’ 
WRITE(*,*)  ’  CGA  =  =  =  >  Enter  0  ’ 

WRITE(*,*)  ’  EGA  ==  =  >  Enter  1  ’ 

READ(*,*)  NS 
WRITE(*,*) 

WRITE{*,*)  ’Enter  number  for  'Form  Feed*  after  plot’ 

WRITE(*,*)’  Enter  0  ===>  NO  FORMFEED  (Laser)’ 

WRITE(*,*r  Enter  1  ===>  FORMFEED  (Impact)’ 

WRITE(*,*j 

READ(^*)  NFD 

ILINE=2 

Pi=3.14159 

C0=2.9979E+8 

Tht=0.0 

Th=0 

a1=(-1.0)*LOG(0.1) 

WRITE(*,*) 

WRITE(*,*)  ’  ENTER  :  String  length  in  Meters’ 
WRiTE(*,*) 

READ(*,*)  L 
WRITE(*,*) 


PROGRAM  TH7.FOR  COHEN  YUVAL  MAY  1990  ’ 
SPACE-TIME  WAVE  EQUATION  PROGRAM 

The  program  solves  numerically  the  displacement’ 
of  a  finite  slring.The  equation  for  U(x,t)  is  the  ’ 
wave  equation  including  a  term  of  loss.  ’ 

The  program  uses  the  "Finite  Difference  Method*  ’ 
to  solve  the  Partial  Differential  Equation  (PDE)  ’ 

The  string  with  length  L  meter  is  suodivided  into  ’ 

N  segments.  The  boundary  conditions  are  ’ 
U(x,t)=U(L,t)=0  and  initial  conditions  U(x,0)=0  ’ 
and  dU(x,o)/clt=0.  The  excitation  is  a  gaussian  ’ 
pulse:  G=Ga*EXP(-a*(t-Tmax-x*tan(Theta))**2)  ’ 

The  pulse  width  is  set  between  *10%*  points.  The  ’ 
incident  angle  is  Theta.  The  space  step  is  S  and  ’ 
the  time  step  is  T.’ 

Ready  to  begin  ? 
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WRITE(*,*) 

WRITE(*,*)  ’  ENTER  :  Gaussian  pulse  parameters’ 
WRITE(*,*)  '  Enter:  Amplitude  GA’ 

WRITE(*,*) 

READ(**)  GA 
WRITER,*) 

WRITEC**) 

WRITER*) 

Tl=1e-9 

T2=0.6*T1 

c 

CONTINUE 

NNE)CT=0 

WRITE(*,*)  ’  ENTER  number  of  segments  N’ 

WRITE(*,*)  ■  (Any  integer  between  2  to  100)’ 

READ(**)  N 
WRITE(*.*) 

S=U/N 

c  *  choose  S/T=C’ 

T=S/C0 

WRITE(*,*)  ’  How  many  TIME  steps  to  compute  ?’ 
WRfTE(*,*)  ’  Enter  Integer  between  2  to  400’ 

WRITE(**) 

READ(**)  M 
WRITE(*,*) 

WR1TE(*,*)  ’  Enter :  coefficient  R  ’ 

WRITE(*,*)  ’  ( R>0  or  R=0 )  in  IE-10  Units’ 

WRrrE(V)  ’  The  case  R=0  is  LOSSLESS  case.’ 
WRITE(*,*) 

READ(*,*)  R 
R=R*1E-10 
Cl  =C0*R»S 
A=2/(2+C1) 

P=(Cl/(2-fCl))-A 

D=A*(S**2)*(-1.0) 

IF(NNEXT.EQ.O)  GO  TO  8 
c 

5  CONTINUE 

WRITER,*)  ’  Enter  Pulse  width  in  NSEC’ 

READ(**)  T1 
Tl=Tl*1E-9 

WRrrE(*,*)  ’  Enter  Time  delay  in  NSEC  for  Peak  Pulse  Impact’ 
READ(*,*)  T2 
T2=T2*1E-9 
NNE)CT=0 
8  CONTINUE 
Tmax=T1/2 
a2=al/CTmax**2) 

IFCT2.GT.Tmax)  GO  TO  9 

WRrrE(*,*)  ’  Time  delay  must  be  longer  then  half  pulse  width!' 
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TRY  AGAIN  !!!’ 


WRITE(*,*)  ’ 

GOTO  5 

6  CONTINUE 

7  CONTINUE 

WRITE(* .*)  ’  Enter  Incident  Angle  in  DEGREE  ( 0<  =Theta=<90 )’ 
READ(**)  Th 

IF  (Th  .GT.  90  .or.  Th  .LT.  0)  then 

WRITE(**)  ’  TRY  AGAIN  !  0=<  Theta  <=90  ’ 

GOTO? 

END  IF 
NNE)CT=0 
9  CONTINUE 
c 

c  GENERATE  BOUNDARY  and  INITIAL  CONDITIONS 
D0  300j=1,M+2 
DO  400  i=1,N+1 
U(i.]-)=0. 

400  CONTINUE 
300  CONTINUE 
c 

c  GENERATE  GAUSSIAN  PULSE  EXCITATION  E=G*EXP(-a*(t-Tmax)**2) 
c  PULSE  WIDTH  at  '1 0%’  POINTS, 
c 

Tht=Crh/180.0)*Pi 

T3=T*tan(tht) 

DO  600  1=1, N 
DO  500  j=1,M+2 

F{i,D=GA*exp((-1.O)*A2*(0-1)*T-T2-((i-1)*T3))**2) 

500  CONTINUE 
600  CONTINUE 

WRITE(*.*)  ’  Want  to  create  file  f.maf  with  driver  data  ?’ 

WRITE{V)’  NO  =  =  >  Enter  0’ 

WRITER, *)’  YES  =  =  >  Enter  1’ 

READC,*)  I 

IF  (I  .EQ.  0)  GO  TO  23 

OPEN(4,file=’f.mat') 

DO  100  j=1,M+2 
DO  200  i=1,N 
WRITE(4.*)  f(i.j) 

200  CONTINUE 
100  CONTINUE 
CLOSE  (4) 

23  CONTINUE 
c 

c  DISPLACEMENT  COMPUTATION  :  MARCHING  FORWARD  IN  TIME 
c 

DO  700  j=3.M+2 

ij=i-2 

DO  800  i=2,N 
ii=i-1 
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Ul=A*{U(i+1,j-1)+U{i-1.j.1)) 

U2=P*U(i,j-2)+D*F(ii.jj) 

U(i,j-)=U1+U2 
800  CONTINUE 
700  CONTINUE 
C 

c  CREATE  a  DATA  FILE  for  OUTPUT  POST  PROCESSING  (matlab) 
WRITE{*,*)  ’  Want  to  create  file  ‘U-mat*  with  displacement  data  ?’ 
WRITE{V)’  NO  ==>  Enter  O’ 

WRITE(**)’  YES  ==>  Enter  1* 

READ(**)  I 

IF  (I  .EQ.  0)  GO  TO  24 

OPEN(3,file=*U.mat’) 

DO  710i=2,M+2 
D0  720i=1,N+1 
WRITE(3 ,*)  UO.D 
720  CONTINUE 
710  CONTINUE 
CLOSE  (3) 
c 

C  OUTPUT 
c 

24  CONTINUE 

IF(NNEXT.EQ.O)  GO  TO  10 

2  CONTINUE 
WRITER, 20)  M 

20  FORMATC  Enter  TIME  in  Time  Step  Units:  It=1,2...M=’,l3,T) 
READ(*  *)  J 

OPEN(2.FILE=’X') 

WRITE(2,11)  J.N.Th 

1 1  FORMAT{'Displacement  at  TIME  t=’,i3,’  STEPS.  ’,i3, 

S’  SEG.,  Theta=’,i3) 

REWIND(2) 

READ(2,12)  TITLE 

12  FORMAT(A) 

REWIND{2) 

NPTS=N+1 

Xmin=0 

Xmax=N 

j=j+2 

DO  900  i=1  ,  N+1 
PL(D=U{i.D 
900  CONTINUE 

CALL  PLTSUB(title,npts.xmin,xmax,xmin,xmax,pI,ns,nfd,iline) 
GOTO  10 
c 

3  CONTINUE 
N1=N-1 
WRITE(*,21)  N1 

21  FORMATC  Enter  SEGMENT  NUMBER:  (1,2, ...N-1  =’.13.’]’) 
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READ(* ,*)  i 
0PEN(2,FILE=’X’) 

WRITE(2,13)  i,N.Th 

13  FORMAT(’DISPLACEMENT  of  SEGMENT  #'.13,'.  M3.’  SEG.,  Theta=’,i3) 
REWIND(2) 

READ(2.14)  TITLE 
REWIND(2) 

14  F0RMAT(A64) 
i=i+1 
NPTS=M 
Xmin=0 
Xmax=M 

DO  910j=1.M+1 
PLOUO.i+1) 

910  CONTINUE 

CALLPLTSUB(titIe,npts,xmin,xmax,xmln,xmax,pl,ns,nfd,iline) 

GOTO  10 

4  CONTINUE 
N1=N-1 

WRrrE(*.22)  N1 

22  FORMATC  Enter  SEGMENT  NUMBER:  [1,2,...N-1  =*.13,’] ’) 

READ(*,*)  i 
OPEN(2,FILE=’X’) 

WRrrE(2,15)  i.T11.T22,N,Th 

1 5  FORMATCDriver  on  Seg.:',i3.’.Tpw=’.F6.3,'  ns.Td=’,F6.3, 

$■  ns,N=M3,’.Theta=’,i3) 

REWIND(2) 

READ(2,16)  TITLE 
REWIND(2) 

16  FORMAT(A64) 

NPTS=M 

Xmin=0 

Xmax=NPTS 

DO920j=1,NPTS 

920  CONTINUE 

CALL  PLTSUB{tftle,npts,xmin,xmax,xmin,xmax,pI,ns,nfd,iTine) 

10  CONTINUE 
WRITER,*) 

WRITE(*,*) '  SELECT  Number  for  Results  or  Change  data' 


WRITE(*,*)’ - ’ 

WRITEC,*)  ’  Change  Data  AGAIN??  ===>  1’ 

WRITE(V)  ’  Displacement  on  the  wire  at  time  t=..  ===>  2' 

WRITE(*,*) '  Displacement  on  Segment  #..  ===>3* 

WRITE{*  *)  ’  Pulse  Excitation  as  Function  of  Time  ===>  4’ 

WRrTE{*,*) '  Change  Pulse  Excitation  Timing  ===>  s' 

VWilTE(V)  ’  Change  Pulse  Excitation  Incident  ANGLE  ===>  6’ 

T11=T1*lE+9 
WRITE{*.*) 
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WRITE(*,17)T11 

17  FORMAT  ('  Currfint '/alues  are:  Pulse  wiclth=’,f6.3,’ ns.’) 
T22=T2*1E+9 

WRITE(*,18)  T22 

18  FORMAT  ( ’  Time  Delav  for  Peak  Pulse  lmpact=’,F6.3,’  ns.’) 
WRiTE(M9)  Th 

19  FORMAT  ( ’  Incidem  Angle  :  Theta  =’,13,’ Degree.’) 

WRITE(^*) 

WRITEC^,*)  ’  Any  Oi:..?r  Integer  ===>  END  PROGRAM  !!!’ 
WRITE(*,*) 

READf,*)  NNC.’^T 
IF{NNEXT.EQ.1)  GO  TO  1 
IF(NNEXT.EQ.2)  GO  TO  2 
IF(NriEXT.EQ.3)  GO  TO  3 
IF(NNE>Cr.EQ.4)  GO  TO  4 
IF(NNEXT.EQ.5)  GO  TO  5 
IF(NNEXT.EQ.6)  GO  TO  6 
STOP 
END 


c 

SUBROUTINE  PLTSUBCTITLE,NPTS,XMI1n!,XMAX,XS,XF,F,NS,NFD,ILINE) 
C 

C  MS-FORTraN  Subroutine  using  "Grafmatic*  Library  Subroutines. 

C  Solid  Line  Using  Portions  of  "PLOP  Program. 

0  Written  by  M.A.  Morgan  with  Latest  Update  to  EGA/CGA  June  1989. 

C 

C  Default  Printer  is  'IBM  Graphics'  (e.g.  Epson,  Okidata,  IBM) 

C  With  Plot  Rotated  90  degrees  From  the  Vertical.  'GrafPlus.Com' 

C  May  be  Run  tc  Rotate  Plot  Upright  on  Paper  and  to  Use  a  Variety 

C  of  Impact  Printers.  'GrafLaser.Com'  May  be  Run  to  Use  a  Laser 

C  Printer.  See  GrafPlus/Laser  Manual  From  Jewell  Technology. 

C 

C  INPU,  DATA  FORMAT 

C 

C  TITLE  -  64  Cha-'acter  Header 

C  MPTS  -  #  Data  Points  id 

C  XMiN  -  Min  X-Axis  value 

C  XMAX  -  Max  X-Axis  value 

C  XS  -Sturti'ig  X  value  for  F(1) 

C  XF  -  FinisI  ling  X  value  for  F(N) 

C  F(N)  -  Input  Data  Array 

C  NS  -  Monitor:  '0'=  CGA,  *1'=  EGA 

C  NFD  -  '1'  ~>  Form  Feed  After  Plot  Hardcopy  (Impact) 

C  Any  Oth.?r  Integer  ->  No  Form  Feed  (Laser) 

C  ILINE  -  'I*  ->  +  +  -f  Symbol  Dot  Plot 

C  Any  Other  Integer  ->  Solid  Line  Plot 

C 
C 
C 


NOTE:  The  X-Axis  Range  Can  Be  Made  Larger  Than 
The  Domain  of  the  Plotted  Function,  F(x). 


o  o 


otherwise,  Make  XS=XMIN  and  XF=XMAX. 

CHARACTER*!  DUM, BELL, FEED 
CHARACTER*64  TITLE,  FN  AM  E,HCOPY 
REAL  X(512),F(512) 

INTEGER*2  N,JROW,JCOL,ISYM,ITYPE,NSCRN 
iNTEGER*2CYAN,WHITE,YELLOW,RED,BLACI<,BLUE,NTWO 
INTEGER*2  JROW1,JROW2,JCOL1,JCOL2 
BELL=CHAR(7) 

FEED=CHAR(12) 

WRITE(*,*)  BELL 
WHITE=7 
CYAN=11 
YELL0W=14 
RED=12 
BLACK=0 
BLUE=1 
NTW0=2 

NSCRN  =  6  +  10*NS 
IF(NFD.EQ.I)  •jPtN(1,FILE=:’PRN’) 

N=NPTS 

DX=(XF-XS)/(N.1.0) 

FMIN=0.0 
FMAX=0.0 
D0  22K=1,N 
X(K)=XS+(K-1.0)*DX 
IF(F(K).LT.FMIN)  FMIN=F(K) 

22  IF(F(K).GT.FMAX)  FMAX=F(K) 

IF(FMIN.GT.O.O)  FMIN=0.0 
IF(FMAX.LT.0.0)  FMAX=0.0 
0  Computing  Scale  Factors  for  Vertical  Axis 
ABSMIN=ABS(FMIN) 

ABSMAX=ABS(FMAX) 

YBIG=AMAX1  (ABSMlN,ABSMAX) 

N3CL=INT{LOG10(YBIG)) 

IF  (YBIG.LT.i.O)  NSCL=NSCL-1 
YSCL=10.**NSCL 
FMlN=FMIN/YSCL 
FMAX=FMAXA'SCL 
ABSMIN=ABSMINA'SCL 
ABSMAX=ABSMAXA'SCL 
D0  33K=1,N 
33  F(K)=F(K)A'SCL 
C  Adaptive  Scaling  Algorithm  5/89 
C  Setting  Polarity  of  YMIN  and  YMAX 
YMIN=0.0 
YBIG=0.0 

IF(FMIN.EQ.O.O)  GO  TO  37 
DY=0.5 
35  CONTINUE 
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IF(YMIfi.GE.4.0)  DY=1.0 
YMIN=YMIN+DY 
IF(ABSMlN.GT.YMIN)  GO  TO  35 
YBIG=YMIN 

YMIN=YMIN*FMIN/ABSMIN 
37  YMAX=0.0 

IF(FMAX.EQ.O.O)  GO  TO  41 
DY=0.5 

39  CONTINUE 

IF(YMAX.GE.4.0)  DY=1.0 
YMAX=YMAX+DY 
!F(ABSMAX.GT.YMAX)  GO  TO  39 
IF(YMAX.GT.YMIN)  YBIG=YMAX 
YMAX=YMAX*FMAX/ABSMAX 
41  CONTINUE 

C  Calling  GRAFMATIC  Routines  and  Plotting  F 
C  Solid  Line  Graph  Default 
ITYPE=1 
ISYM=-2 
NDOTS=0 

C  +  +  +  Line  Graph  If  ILINE=1 
IF(ILINE.NE.1)  GO  TO  45 
ISYM=-1 
ITYPE=0 
45  CONTINUE 

JROW1=  35+25*NS 
JROW2=  170+120*NS 
JCOL1=  110-10*NS 
JCOL2=  640 
CALL  QSMODE(NSCRN) 

CALL  QPTXT(64, TITLE, YELLOW, 16,24) 

CALL  QPLOT{JCOL1 , JCOL2, JROW1  ,JROW2,XMIN,XMAX,YMIN,YMAX,XMIN, 
1  0..0,1.,1.5) 

CALL  QSETUP(NDOTS,CYAN,ISYM.CYAN) 

XMAJOR= (XM  AX-XMIN)/5.0 

CALL  QXAXIS(XMIN,XMAX,XMAJOR,1,1,2) 

IF(YBIG.LE.4.0)  YMAJOR=0.5 
IF{YBIG.GE.5.0)  YMAJOR=1.0 
IF(YBIG.EQ.8.0)  YMAJOR=2.0 
IF(YBIG.EQ.9.0)  YMAJOR=:3.0 
IF(YBIG.EQ.10.)  YMAJOR=2.0 
CALL  QYAXIS(YMIN,YMAX,YMAJOR.1,1,1) 

JROW=32+21*NS+(ABS(YMIN)/(ABS(YMAX)+ABS(YMIN)))*(135+95*NS) 

JCOL=80-8*NS 

CALL  QGTXT(3,’0.0’,WHITE,JCOL,JROW,0) 

CALL  QPTXT(1,’S’,YELLOW,5,18) 

CALL  QPTXT(1,’c’,YELLOW,5,17) 

CALL  QPTXT(1,’a’,YELLOW,5,16) 

CALL  QPTXT(1,T,YELLOW,5,15) 

CALL  QPTXT(1,’e’, YELLOW, 5,14) 
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CALL  QPTXT(1,’X’,YELLOW.5.12) 

CALL  QPTXT(2,’1  O’, YELLOW, 4,10) 

CALL  QPTXT(1,’*’,YELLOW,5,9) 

CALL  QPTXT(1,’*’,YELLOW,5,8) 

CALL  QCMOV(0,8) 

WRITE(*,150)  NSCL 
CALL  QTABL(ITYPE,N,X.F) 

C  Hardcopy  Query  -  Updated  5/1 1/89 
HCOPY=’Hardcopy  —>  Enter  P  or  p’ 

CALL  QPTXT(30,HCOPY.RED,25,1) 

CALL  QCMOV(55,1) 

READ(*,100)  DUM 
HCOPY=’ 

CALL  QPTXr(40,HCOPY.BLACK,25,1) 
iF(DUM.NE.’P’  .AND.  DUM.NE.’p’)  GO  TO  48 
CALL  QPSCRN 

IF(NFD.EO.I)  WRITE(1,160)  FEED 
48  CONTINUE 

C  Exit  to  Blue  Background  on  Screen  -  To  Change  This, 

C  Replace  'BLUE*  in  Calls  to  QPREG  and  QOVSCN  to  That  Desired. 
CALL  QSMODE(NTWO) 

CALL  QPREG(0,BLUE) 

CALL  QOVSCN(BLUE) 

100  FORMAT(A) 

120  F0RMAT(I5) 

130  FORMAT(E12.3) 

150  F0RMAT(4X,I3.\) 

160  FORMATC  ’,A,\) 

RETURN 

END 
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APPENDIX  B.  ARMA  MODEL  ALGORITHM 

This  Appendix  presents  the  development  of  the  ARMA  model  using  the 
equation  of  evolution  for  the  lossless  case  example  given  in  Chapter  IV. 

Start  with  the  equation  of  evolution  of  the  following  form  for  M=3 

U(k)=A-U(k-l)-U(k-2),  (1) 

with  the  symmetric  matrix 

b  1  0 
4=1  0  1  . 
p  1  p 

Using  Eq.  (1)  to  express  U(k-l),  substituting  into  Eq.  (1)  yields  the  following 

Uik)  =  -  Uik-2)  +A^-U(k-2)  -A'V{k-3) .  (2) 

Using  Eq.  (1)  to  express  U(k-2)  and  substituting  into  Eq.  (2)  yields  the  following 

U{k)  =  -  Uik-2)  +A^'\A'U(k-3) -  Uik-4)] -A-Uik-3) .  (3) 

Rearranging  Eq.  (3)  and  using  the  fact  that  for  M=3  A''=2A,  the  following 
expression  is  obtained: 

Uik)  =  -  Uik-2)  +A-Uik-3)  -^Uik-A) . 
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Using  Eq.  (1)  to  write  U(k-3),  substituting  in  Eq.  (4),  and  adding  and  subtracting 
U(k-4)  yields 


U(k)  =  -  U(k-2) -A'U(k-5)  +  U(k-4)  -  Uik-4) .  (5) 


Using  Eq.  (1)  to  write  U(k-4)  and  substituting  in  Eq.  (5)  yields 


U(k)  =  -  U(k-2)  -  U(k-4)  -  Uik-6) . 


(6) 


Equation  (6)  has  the  form  of  the  ARMA  model  for  M=3  described  in  Chapter  IV 
Section  G. 
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APPENDIX  C.  PEONY’S  METHOD  PROGRAM 


The  computer  program  entitled  TEST.M  implements  the  Prony’s  method  to 
estimate  ARMA  model  coefficients  using  time-domain  data.  The  program,  which  is 
written  using  MATLAB  codes,  finds  the  coefficients  by  implementing  the  procedure 
described  in  Chapter  IV,  Section  1.  The  matrix  with  the  time-domain  data  should  be 
defined  in  a  MATLAB  format,  and  stored  in  a  file  called  X.MAT.  The  data  may  be 
generated  by  TH7.FOR,  and  translated  into  MATLAB  format  using  "translate"  in 
MATLAB.  The  output  of  the  program  plots  the  results  of  the  coefficients  for  each 
segment  data.  The  coefficients  are  stored  in  a  matrix  called  c  where  each  column 
contains  the  coefficients  for  each  segment  data.  Plots  of  results  may  be  obtained  by 
using  the  information  in  mat^’x  c. 


%  TEST.M  program  by  Yuval  Cohen  June  1990 
%  The  program  implements  the  Prony’s  method  on  time-domain  data 
%  estimating  ARMA  coefficients  of  the  late-time  solution  of  the 
%  space-time  wave  equation. 

% 

load  X 

input(’Enter  number  of  segments  of  the  string  N’) 
ns=ans; 

input(’Enter  number  of  coefficients  M=2(N-1)') 
m=ans; 

input(’Enter  the  time-step  f) 
t=ans: 

c=[]: 

for  1=1  :ns  n=0: 

y=[l: 
u=I]: 
for  i=1:m 
for  j=1:m 
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y(i.i)=x(t-i-n,i): 

end 

n=n+1; 

end 

for  i=1:m  u(0=x(t-i+1,l); 

end 

u=u’: 

a=y\u: 

c(:,l)=a: 

plot(a) 

end 
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